*******************************************************************************
*   Title: draft_final.do
* Purpose: Figures and tables for final version of draft 
*   Owner: EZ
*    Date: 2022-07-07
*******************************************************************************

global gpr = "plotregion(fcolor(white) lcolor(white) m(tiny)) graphregion(fcolor(white) lcolor(white))"

global f1 "0 80 164"
global f2 "white"
global f3 "239 65 53"
global f4 "200 200 200"

global u1 "60 59 110"
global u2 "white"
global u3 "178 34 52"
global u4 "50 50 50"
global u5 "205 207 112"
global u6 "100 100 100"

*u1, u3, u5, and u6 are replaced by u7-u10 respectively for select QJE figs
global u7 "60 59 110"
global u8 gs13
global u9 "205 207 112"
global u10 "100 100 100"

global p1 "164 0 162"
global p2 "106 105 195"
global p3 "191 184 246"

*p21 is used to replace p2 for select QJE figs
global p21 "106 105 195"

capture program drop graph_export/*%<*/
program define graph_export
    graph export `1'.pdf, replace
    graph export `1'.eps, replace

end/*%>*/

*******************************************************************************
* Introduction %<
*******************************************************************************
capture program drop graph_comparison_overtime/*%<*/
program define graph_comparison_overtime

    /******************************************************************************
        Wealth Concentration in the United States
    ******************************************************************************/

    /***************************************************************************
        Top 0.1% Share of Total Wealth, Raw Series
    ***************************************************************************/    

    /**********************************************************************
        Prep PSZ equal split numbers, pulling from literature  
            folder in PSZ2017AppendixTablesII(Distrib).xlsx sheet TE2c. 
            Drop observations if they're missing top 0.1% share, rescale  
            into percentages from proportions and save as tempfile.
    **********************************************************************/

    load_analysis_data psz_top
    tempfile pszes
    save `pszes'
                    
    /**********************************************************************
        Prep baseline SZ wealth shares and Kopzcuk-Saez wealth 
            shares, pulling both from SZ QJE 2016 Appendix, in
            AppendixTables(OtherEstimates).xlsx, sheet C4. Save as 
            tempfile.
    **********************************************************************/

    load_analysis_data ks01
    tempfile sz16_ks
    save `sz16_ks'

    /**********************************************************************
        Prep SCF series
    **********************************************************************/

    load_analysis_data scfraw
    reshape wide networth, i(year) j(top01)
    gen total_wealth = networth0 + networth1
    gen top01_scfraw = (networth1 / total_wealth) * 100
    tempfile scfraw
    save `scfraw'

    /***********************************************************************
        Prepare capitalized series
    ***********************************************************************/

    foreach capspec in equrtrns pref {
        
        local hwealname = cond("`capspec'" == "pref", "hweal_preferred", "hweal20")

        if "`capspec'" == "pref" {
            load_analysis_data szz
        }
        else {
            load_analysis_data sz
        }

        keep if inlist(group, "All", "P99.9-100")
        
        sort year group
        by year: assert _N == 2

        replace group = cond(group == "All", "_all", "_top01")

        keep year group `hwealname'

        reshape wide `hwealname', i(year) j(group, string)

        gen top01_`capspec' = (`hwealname'_top01 / `hwealname'_all) * 100

        assert !missing(top01_`capspec')
        assert inrange(year, 1966, 2016)

        tempfile `capspec'
        save ``capspec''
    }

    /**********************************************************************
        Merge together data pulls from different sources by year,  
            and export merged 
    **********************************************************************/

    use `pref', clear
    keep year top01_pref

    merge 1:1 year using `equrtrns', keepusing(top01_equrtrns) assert(3) nogen

    merge 1:1 year using `scfraw', keepusing(top01_scfraw)
    assert year == 2019 if _merge == 2
    assert !inrange(year, 1989, 2019) | mod(year - 1989, 3) > 0 if _merge == 1
    assert _merge == 3 if inrange(year, 1989, 2016) & mod(year - 1989, 3) == 0 
    drop _merge

    merge 1:1 year using `sz16_ks', keepusing(top01_ks top01_sz16)
    assert year > 2012 if _merge == 1
    assert year < 1966 if _merge == 2
    assert _merge == 3 if inrange(year, 1966, 2012)
    drop _merge

    merge 1:1 year using `pszes', assert(1 3) keepusing(top01_pszes)
    assert !missing(year) & year > 2014 if _merge == 1
    assert _merge == 3 if year <= 2014
    drop _merge

    /* Ensure that our extension of PSZ 2018 is actually matching their appendix 
        data, subject to a reasonable imprecision error */
    gen psz18diff = abs(top01_equrtrns - top01_pszes)
    assert psz18diff < 0.42 if !missing(top01_equrtrns) & !missing(top01_pszes)
    drop psz18diff
    
    export delimited using top01shares_unadjusted.csv, replace
            
    /***********************************************************************
        Make graph (2 versions, one for slides)
    ***********************************************************************/

    #delimit ;
    twoway (scatter top01_sz year, c(l) ms(dh) lc("$f3") mc("$f3") lpattern(shortdash)) 
            (scatter top01_pszes year, c(l) ms(s) lc("$u3") mc("$u3")) 
            (scatter top01_equrtrns year, c(l) ms(sh) lc("$u3") mc("$u3") lp(shortdash)) 
            (scatter top01_ks year, c(l) ms(d) lc("$f1") mc("$f1")) 
            (scatter top01_scf year, c(l) ms(T) lc("$p2") mc("$p2") lw(thin))  
            (scatter top01_pref year, c(l) ms(O) lc("$u1") mc("$u1") lwidth(thick)),
            legend(col(2) region(lcolor(white) margin(tiny)) 
                order(1 "Equal Return, Tax Unit (Saez Zucman 2016)" 
                      5 "Raw SCF" 
                      2 "Equal Return, Individual (PSZ 2018)" 
                      4 "Unadjusted Estate Multiplier (Kopczuk Saez 2004)" 
                      3 "Equal Return, Individual (PSZ 2018, Extended)" 
                      6 "Our Baseline Estimate")) 
            $gpr xlab(1915(10)2015) xtitle(" ") 
            ytitle("Share of Total Household Wealth (%)") xsize(8);
    graph_export top01shares_unadjusted

    twoway (scatter top01_sz year, c(l) ms(sh) lc("$f3") mc("$f3") lpattern(shortdash)) 
            (scatter top01_pszes year, c(l) ms(s) lc("$u3") mc("$u3")) 
            (scatter top01_equrtrns year, c(l) ms(sh) lc("$u3") mc("$u3") lp(shortdash)) 
            (scatter top01_ks year, c(l) ms(d) lc("$f1") mc("$f1")) 
            (scatter top01_scf year, c(l) ms(T) lc("$p2") mc("$p2")),  
            legend(col(2) region(lcolor(white) margin(tiny)) 
                order(1 "Equal Return, Tax Unit (Saez Zucman 2016)" 
                      5 "Raw SCF" 
                      2 "Equal Return, Individual (PSZ 2018)" 
                      4 "Unadjusted Estate Multiplier (Kopczuk Saez 2004)" 
                      3 "Equal Return, Individual (PSZ 2018, Extended)" 
                      )) 
            $gpr xlab(1915(10)2015) xtitle(" ") 
            ytitle("Share of Total Household Wealth (%)") xsize(8);
    graph export top01shares_unadjusted_noszz.pdf, replace;
    #delimit cr

    /***************************************************************************
        Panel B: Top 1%, 0.1%, and 0.01% Share of Total Wealth, Harmonized
            Series (feat. preferred, equal return equal split, and harmonized
            SCF).
    ***************************************************************************/    

    /***********************************************************************
        Prepare PSZ 2018 equal-split series
    ***********************************************************************/  

    load_analysis_data psz_top
    tempfile psz18es
    save `psz18es'

    load_analysis_data sz2020_top
    tempfile sz20
    drop if year < 1966
    save `sz20'

    /***********************************************************************
        Prepare post-2014 extension of of equal-returns top wealth 
            shares which we produce on the inside
    ***********************************************************************/  
    load_analysis_data sz

    keep if w${sz_wlth_defn}_group <= 6

    assert inlist(group, "All", "P90-100", "P99-100", "P99.9-100", "P99.99-100", "P99.999-100") 
    isid year group
    
    replace group = cond(group == "P99.999-100", "_top0001", ///
                    cond(group == "P99.99-100", "_top001", ///
                    cond(group == "P99.9-100", "_top01", ///
                    cond(group == "P99-100", "_top1", ///
                    cond(group == "P90-100", "_top10", "_total")))))

    keep year group hweal$sz_wlth_defn
    rename hweal$sz_wlth_defn hweal

    reshape wide hweal, i(year) j(group, string)
    isid year

    foreach topgrp in top10 top1 top01 top001 {
        gen `topgrp'_equrtrns = (hweal_`topgrp' / hweal_total) * 100
        assert inrange(`topgrp'_equrtrns, 0, 100)
    }

    tempfile equrtrns
    save `equrtrns'

    load_analysis_data sz_v3

    keep if w${sz_wlth_defn_late_v3}_group <= 6 | w${sz_wlth_defn_mid_v3}_group <= 6 | ///
        w${sz_wlth_defn_early_v3}_group <= 6

    assert inlist(group, "All", "P90-100", "P99-100", "P99.9-100", "P99.99-100", "P99.999-100") 
    isid year group
    
    replace group = cond(group == "P99.999-100", "_top0001", ///
                    cond(group == "P99.99-100", "_top001", ///
                    cond(group == "P99.9-100", "_top01", ///
                    cond(group == "P99-100", "_top1", ///
                    cond(group == "P90-100", "_top10", "_total")))))

    keep year group hweal_szv3
    rename hweal_szv3 hweal

    reshape wide hweal, i(year) j(group, string)
    isid year

    foreach topgrp in top10 top1 top01 top001 top0001 {
        gen `topgrp'_equrtrns_v3 = (hweal_`topgrp' / hweal_total) * 100
        assert inrange(`topgrp'_equrtrns_v3, 0, 100)
    }

    tempfile equrtrns_v3
    save `equrtrns_v3'

    /***********************************************************************
        Prepare adjusted SCF top shares
    ***********************************************************************/  

    load_analysis_data scfharmony_v3

    tempfile scfpref
    save `scfpref'

    load_analysis_data scfharmony_v3_forbes

    tempfile scfprefforbes
    save `scfprefforbes'

    /***********************************************************************
        Prepare top shares under preferred capitalization
    ***********************************************************************/  
    load_analysis_data szz

    keep if w${preferred_defn_early}_group <= 6 | w${preferred_defn_mid}_group <= 6 | ///
            w${preferred_defn_late}_group <= 6
    
    assert inlist(group, "All", "P90-100", "P99-100", "P99.9-100", "P99.99-100", "P99.999-100") 
    isid year group
    
    replace group = cond(group == "P99.999-100", "_top0001", ///
                    cond(group == "P99.99-100", "_top001", ///
                    cond(group == "P99.9-100", "_top01", ///
                    cond(group == "P99-100", "_top1", ///
                    cond(group == "P90-100", "_top10", "_total")))))
    #delimit cr

    keep year group hweal_preferred
    reshape wide hweal_preferred, i(year) j(group, string)

    foreach topgrp in top10 top1 top01 top001 top0001 {
        gen `topgrp'_preferred = (hweal_preferred_`topgrp' / hweal_preferred_total) * 100
        assert inrange(`topgrp'_preferred, 0, 100)
    }

    /***********************************************************************
        Merge series together
    ***********************************************************************/

    merge 1:1 year using `equrtrns', assert(3) nogen
    merge 1:1 year using `equrtrns_v3', assert(3) nogen

    merge 1:1 year using `scfpref'
    assert year == 2019 if _merge == 2 
    assert !inrange(year, 1989, 2019) | mod(year - 1989, 3) > 0 if _merge == 1
    assert _merge == 3 if inrange(year, 1989, 2016) & mod(year - 1989, 3) == 0
    drop _merge 

    merge 1:1 year using `scfprefforbes', nogen

    merge 1:1 year using `psz18es', keepusing(top*_pszes)
    assert !missing(year) & year > 2014 if _merge == 1
    assert year < 1966 if _merge == 2
    assert _merge == 3 if inrange(year, 1966, 2014)
    keep if inlist(_merge, 1, 3)
    drop _merge

    merge 1:1 year using `sz20', keepusing(top*_pszes20)

    foreach topgrp in top001 top01 top1 top10 { 
        di "`topgrp'"
        gen psz18`topgrp'diff = abs(`topgrp'_equrtrns - `topgrp'_pszes)
        assert psz18`topgrp'diff < 1.01 if !missing(`topgrp'_equrtrns) & !missing(`topgrp'_pszes)
        drop psz18`topgrp'diff
    }

    export delimited using topshares_adjusted.csv, replace

    sort year
    foreach topgrp in top0001 top001 top01 top1 top10 {

        if ("`topgrp'" == "top0001") {
            local yscale "yscale(range(0.5 5))"
        }
        else if ("`topgrp'" == "top001") {
            local yscale "yscale(range(1.5 10))"
        }
        else if ("`topgrp'" == "top01") {
            local yscale "yscale(range(5.5 20))"
        }
        else if ("`topgrp'" == "top1") {
            local yscale "yscale(range(20.5 40))"
        }
        else if ("`topgrp'" == "top10") {
            local yscale "yscale(range(58 78))"
        }
        
        #delimit ;
        graph twoway 
            (connect `topgrp'_equrtrns year, ms(sh) lc("$u3") mc("$u3") lpattern(-)) 
            (connect `topgrp'_pszes20 year, ms(s) lc("$u3") mc("$u3"))
            (connect `topgrp'_pref year, ms(O) lc("$u1") mc("$u1") lwidth(thick))
            (connect `topgrp'_scfpref year, ms(T) lc("$p2") mc("$p2") lwidth(thin))
            (connect `topgrp'_scfpref_forbes year, ms(Th) lc("$p2") mc("$p2") lwidth(thin))
            ,
            $gpr
            legend(order(
                  1 "Equal Return, Individual (PSZ 2018, Extended)" 
                  2 "Revised Saez Zucman (2020)"
                  3 "Our Baseline Estimate"
                  4 "Harmonized SCF" 
                  5 "Harmonized SCF w/Forbes" )
                row(5) region(lcolor(white) margin(tiny)))
            xlab(1960(10)2020) xtitle(" ")
            `yscale'
            ytitle("Share of Total Household Wealth (%)") xsize(4);
        #delimit cr
        graph_export "`topgrp'shares_adjusted"
    }

end/*%>*/
* %>
*******************************************************************************

*******************************************************************************
* Hetero returns%<
*******************************************************************************
****************************************************************************
* Fixed 
****************************************************************************
capture program drop graph_scf_portfolio_fixed/*%<*/
program define graph_scf_portfolio_fixed

    /***************************************************************************
        SCF fixed income portfolio shares
    ***************************************************************************/
      
    /**********************************************************************
        Load data and rank by preferred wealth group; then 
            partition into wealth groups
    **********************************************************************/
    load_analysis_data scfcross_v3

    keep if year == 2016

    * Group has cuts at 0, .9(.01).99 .999 .9999 1.01
    gen wlthgrp = cond(grp > 10, 4, ///
                  cond(grp > 9, 3, ///
                  cond(grp > 0, 2, 1)))

    /**********************************************************************
        Interested in portfolios of fixed claims excluding money 
            market funds and munis.
    **********************************************************************/
    gen hwfix_asset_excl_muni = hwfix_scf //- intexmw

    gen currency_and_deposits = mmda + saving + call + cds + currency + mmmf

    #delimit ;
    gen bonds_fixmutf_excl_mmmf = 
              /* Taxable Bonds */ savbnd + (bond - notxbnd) + 
               /* Mutf excl MM */ (0.5 * comutf) + gbmutf + obmutf + (0.5 * omutf) + 
                      /* Loans */ privloans + mortgageassets + 
      /* Trusts in txble fixed */ trusts_inttaxw + trusts_mmbondfund + 
           /* Tax-exempt bonds */ intexmw;                   
    #delimit cr

    gen checktotal = currency_and_deposits + bonds_fixmutf_excl_mmmf

    assert inrange(hwfix_asset_excl_muni / checktotal, 0.999, 1.001) | ///
           hwfix_asset_excl_muni == checktotal
    drop checktotal

    /**********************************************************************
       Collapse to yield totals by wealth group and calculate
            wealth group specific shares
    **********************************************************************/
    collapse (sum) hwfix_asset_excl_muni currency_and_deposits ///
        bonds_fixmutf_excl_mmmf, by(wlthgrp)

    foreach hwfixcmpt of varlist currency_and_deposits bonds_fixmutf_excl_mmmf {
        gen `hwfixcmpt'_sh = (`hwfixcmpt' / hwfix_asset_excl_muni) * 100
    }

    egen total_check = rowtotal(*_sh)
    assert abs(total_check - 100) < 0.01
    drop total_check

    /**********************************************************************
        Make graph(s)
    **********************************************************************/

    #delimit; 
	
	
    graph bar (asis) currency_and_deposits_sh bonds_fixmutf_excl_mmmf_sh
        ,
        bargap(20) bar(1, color("$u7")) bar(2, fcolor("$u8") lcolor(black))
        legend(order(1 "Currency and deposit-type investments" 
                     2 "Bonds and fixed income funds (excl. money market)")
               region(lcolor(white) margin(tiny)) row(2))  
        over(wlthgrp, relabel(1 "P0-90" 2 "P90-99" 3 "P99-99.9" 4 "P99.9-100")) 
        $gpr
        ytitle("Share of Fixed Income Claims (%)") xsize(5.5);
    #delimit cr
    graph_export "hwfix_asset_excl_muni_portfolio_heterogeneity"
end/*%>*/

capture program drop graph_interest_participation/*%<*/
program define graph_interest_participation

    load_analysis_data interest_info
    keep if year == 2016

    isid wadjgross

    foreach intsrc in "1099bank" "1099sav" "1099loan" "1065" "1120s" "1041" {
        gen share_`intsrc' = 100 * (has_`intsrc' / count)
    }

    #delimit ;
    twoway 
        (connect share_1099bank wadjgross, msize(small) mc("$u1") lc("$u1") lw(thin) ms(o)) 
        (connect share_1099loan wadjgross, msize(small) mc("$p2") lc("$p2") lw(thin) ms(d))
        (connect share_1120s wadjgross, msize(small) mc("$f3") lc("$f3") lw(thin) ms(t)) 
        (connect share_1041 wadjgross, msize(small) mc("$u4") lc("$u4") lw(thin) ms(+) lp("-"))
        (connect share_1099sav wadjgross, msize(small) mc("$f1") lc("$f1") lw(thin) ms(oh))
        (connect share_1065 wadjgross, msize(small) mc("$u3") lc("$u3") lw(thin) ms(s)) 
        ,
        xtitle("AGI Percentile") ytitle("Share with Positive Income (%)")
        legend(order(1 "1099-INT Bank" 6 "1065-K1" 3 "1120S-K1" 
                     5 "Savings Bond" 2 "1099-INT Loan" 4 "1041-K1") col(3)
               region(lc(white))) $gpr xsize(5.5);
    #delimit cr
    graph_export "interest_participation_adjgross_2016"

end/*%>*/

capture program drop graph_interest_composition/*%<*/
program define graph_interest_composition

    /***************************************************************************
       Sources of taxable interest income by AGI rank (lifted straight 
            from Eric's code in xpez_20200117.do, in the repo)
    ***************************************************************************/

    load_analysis_data interest_info
    keep if year == 2016

    drop *_sec *_kid*
    collapse (sum) interest_* txintamt adjgross, by(wadjgross)

    foreach intsrc in "1099bank" "1099sav" "1099loan" "1065" "1120s" "1041" {
        gen incomecompo_`intsrc' = 100 * (interest_`intsrc' / txintamt)
    }

    egen incomecompo_info = rsum(incomecompo_1099bank-incomecompo_1041)

    foreach intsrc in "1099bank" "1099sav" "1099loan" "1065" "1120s" "1041" {
        gen incomecompo_`intsrc'R = 100 * incomecompo_`intsrc'/incomecompo_info
    }

    gen incomecompo_na = 100 - incomecompo_1099bank - incomecompo_1099loan - ///
                         incomecompo_1099sav - incomecompo_1065 - incomecompo_1120s - ///
                         incomecompo_1041

    /***********************************************************************
        Top decile of AGI distribution only
    ***********************************************************************/

    #delimit ;
    twoway 
        (connect incomecompo_1099loanR wadjgross, msize(small) mc("$p2") lc("$p2") lw(thin) ms(d))
        (connect incomecompo_1120sR wadjgross, msize(small) mc("$f3") lc("$f3") lw(thin) ms(t)) 
        (connect incomecompo_1041R wadjgross, msize(small) mc("$u4") lc("$u4") lw(thin) ms(+) lp("-"))
        (connect incomecompo_1099savR wadjgross, msize(small) mc("$f1") lc("$f1") lw(thin) ms(oh))
        (connect incomecompo_1065R wadjgross, msize(small) mc("$u3") lc("$u3") lw(thin) ms(s)) 
        (connect incomecompo_1099bankR wadjgross, msize(small) mc("$u1") lc("$u1") lw(thin) ms(o)) 
            if wadjgross >= 75,
        ysca(range(-2.5 65)) ylab(0(10)60)
        xsca(range(74.5 100.5)) xlab(75(5)100)
        xtitle("AGI Percentile") ytitle("Share of Bin's Taxable Interest (%)")
        legend(order(6 "1099-INT Bank" 5 "1065-K1" 2 "1120S-K1" 
                     4 "Savings Bond" 1 "1099-INT Loan" 3 "1041-K1") col(3)
               region(lc(white))) $gpr xsize(5.5);
    #delimit cr
    graph_export "interest_incomecompo_p90_adjgross_2016"

end/*%>*/

capture program drop graph_bankshares_overtime/*%<*/
program define graph_bankshares_overtime

    load_analysis_data interest_info
    gen share_pos1099bank = (has_1099bank / count) * 100
    gen share_posnonqual = (has_nonqualdivs / count) * 100

    keep year wadjgross share_pos1099bank share_posnonqual
    reshape wide share_pos1099bank share_posnonqual, i(wadjgross) j(year)

    sort wadjgross
    colorpalette "$u1", intensity(0.1(.05)1)
    #delimit ;
    twoway 
        (connect share_pos1099bank2004 wadjgross, ms(o) msize(small) mc("`r(p2)'") lc("`r(p2)'")) 
        (connect share_pos1099bank2006 wadjgross, ms(o) msize(small) mc("`r(p5)'") lc("`r(p5)'"))
        (connect share_pos1099bank2008 wadjgross, ms(o) msize(small) mc("`r(p8)'") lc("`r(p8)'"))
        (connect share_pos1099bank2010 wadjgross, ms(o) msize(small) mc("`r(p11)'") lc("`r(p11)'")) 
        (connect share_pos1099bank2012 wadjgross, ms(o) msize(small) mc("`r(p14)'") lc("`r(p14)'")) 
        (connect share_pos1099bank2014 wadjgross, ms(o) msize(small) mc("`r(p17)'") lc("`r(p17)'")) 
        (connect share_pos1099bank2016 wadjgross, ms(o) msize(small) mc("`r(p19)'") lc("`r(p19)'")) 
        ,
        $gpr
        xtitle("AGI Percentile") ytitle("Share with 1099-INT Bank Income (%)")
        yscale(range(-0.5 100.5))
        legend(order(1 "2004" 2 "2006" 3 "2008"
                     4 "2010" 5 "2012" 6 "2014" 7 "2016") row(2)
               region(lc(white)) margin(tiny)) xsize(5.5);
    #delimit cr
    graph_export "particip_1099bank_yeargradient"

    colorpalette "$u3", intensity(0.1(.05)1)
    #delimit ;
    twoway  
        (connect share_posnonqual2004 wadjgross, ms(s) msize(small) mc("`r(p2)'") lc("`r(p2)'")) 
        (connect share_posnonqual2006 wadjgross, ms(s) msize(small) mc("`r(p5)'") lc("`r(p5)'"))
        (connect share_posnonqual2008 wadjgross, ms(s) msize(small) mc("`r(p8)'") lc("`r(p8)'"))
        (connect share_posnonqual2010 wadjgross, ms(s) msize(small) mc("`r(p11)'") lc("`r(p11)'")) 
        (connect share_posnonqual2012 wadjgross, ms(s) msize(small) mc("`r(p14)'") lc("`r(p14)'")) 
        (connect share_posnonqual2014 wadjgross, ms(s) msize(small) mc("`r(p17)'") lc("`r(p17)'")) 
        (connect share_posnonqual2016 wadjgross, ms(s) msize(small) mc("`r(p19)'") lc("`r(p19)'")) 
        ,
        $gpr
        xtitle("AGI Percentile") ytitle("Share with Non-Qualified Dividends (%)")
        yscale(range(-0.5 100.5))
        legend(order(1 "2004" 2 "2006" 3 "2008"
                     4 "2010" 5 "2012" 6 "2014" 7 "2016") row(2)
               region(lc(white)) margin(tiny)) xsize(5.5);
    #delimit cr
    graph_export "particip_nonqualdivs_yeargradient"

    colorpalette "$u1", intensity(0.1(.05)1)
    forv i = 1/19 {
        local p`i' = "`r(p`i')'"
    }

    colorpalette "$u3", intensity(0.1(.05)1)
    #delimit ;
    twoway  
        (connect share_posnonqual2004 wadjgross, ms(s) msize(small) mc("`r(p2)'") lc("`r(p2)'")) 
        (connect share_posnonqual2005 wadjgross, ms(s) msize(small) mc("`r(p3)'") lc("`r(p3)'")) 
        (connect share_posnonqual2006 wadjgross, ms(s) msize(small) mc("`r(p5)'") lc("`r(p5)'"))
        (connect share_posnonqual2007 wadjgross, ms(s) msize(small) mc("`r(p7)'") lc("`r(p7)'")) 
        (connect share_posnonqual2008 wadjgross, ms(s) msize(small) mc("`r(p8)'") lc("`r(p8)'"))
        (connect share_posnonqual2009 wadjgross, ms(s) msize(small) mc("`r(p10)'") lc("`r(p10)'")) 
        (connect share_posnonqual2010 wadjgross, ms(s) msize(small) mc("`r(p11)'") lc("`r(p11)'")) 
        (connect share_posnonqual2011 wadjgross, ms(s) msize(small) mc("`r(p13)'") lc("`r(p13)'")) 
        (connect share_posnonqual2012 wadjgross, ms(s) msize(small) mc("`r(p14)'") lc("`r(p14)'")) 
        (connect share_posnonqual2013 wadjgross, ms(s) msize(small) mc("`r(p16)'") lc("`r(p16)'")) 
        (connect share_posnonqual2014 wadjgross, ms(s) msize(small) mc("`r(p17)'") lc("`r(p17)'")) 
        (connect share_posnonqual2015 wadjgross, ms(s) msize(small) mc("`r(p18)'") lc("`r(p18)'")) 
        (connect share_posnonqual2016 wadjgross, ms(s) msize(small) mc("`r(p19)'") lc("`r(p19)'")) 
        (connect share_pos1099bank2002 wadjgross, ms(t) msize(small) mc("`p1'") lc("`p1'")) 
        (connect share_pos1099bank2003 wadjgross, ms(t) msize(small) mc("`p2'") lc("`p2'")) 
        (connect share_pos1099bank2004 wadjgross, ms(t) msize(small) mc("`p3'") lc("`p3'")) 
        (connect share_pos1099bank2005 wadjgross, ms(t) msize(small) mc("`p4'") lc("`p4'")) 
        (connect share_pos1099bank2006 wadjgross, ms(t) msize(small) mc("`p5'") lc("`p5'"))
        (connect share_pos1099bank2007 wadjgross, ms(t) msize(small) mc("`p6'") lc("`p6'"))
        (connect share_pos1099bank2008 wadjgross, ms(t) msize(small) mc("`p7'") lc("`p7'"))
        (connect share_pos1099bank2009 wadjgross, ms(t) msize(small) mc("`p9'") lc("`p9'"))
        (connect share_pos1099bank2010 wadjgross, ms(t) msize(small) mc("`p11'") lc("`p11'")) 
        (connect share_pos1099bank2011 wadjgross, ms(t) msize(small) mc("`p13'") lc("`p13'"))
        (connect share_pos1099bank2012 wadjgross, ms(t) msize(small) mc("`p14'") lc("`p14'")) 
        (connect share_pos1099bank2013 wadjgross, ms(t) msize(small) mc("`p16'") lc("`p16'"))
        (connect share_pos1099bank2014 wadjgross, ms(t) msize(small) mc("`p17'") lc("`p17'")) 
        (connect share_pos1099bank2015 wadjgross, ms(t) msize(small) mc("`p18'") lc("`p18'"))
        (connect share_pos1099bank2016 wadjgross, ms(t) msize(small) mc("`p19'") lc("`p19'")) 
        ,
        $gpr
        xtitle("AGI Percentile") ytitle("Share with Positive Income (%)")
        yscale(range(-0.5 100.5))
        legend(order(28 "1099-INT Bank (2002-2016)" 13 "Non-Qualified Dividends (2004-2016)") row(1)
               region(lc(white)) margin(tiny) symxsize(*.5)) xsize(5.5);
    #delimit cr
    graph_export "particip_banknqd_yeargradient"

end/*%>*/

capture program drop graph_asset_returns/*%<*/
program define graph_asset_returns

    /***************************************************************************
    ***  Asset Class Returns [PLACEHOLDER PENDING MORE EZ WORK] ***********
    ***************************************************************************/

    /***********************************************************************
        Retrieve average Savov deposit rate for 2016 and bank rates
        implied by info returns capitalization
    ***********************************************************************/
    load_analysis_data deposit_rates
    sum r_deposit if year == 2016
    local bankrate = `r(mean)'

    load_analysis_data bank_rates
    tempfile bank
    save `bank'

    /***********************************************************************
        Retrieve different estimates of boutique interest rates 
            along the AGI distribution from EZ work, and saving bonds and loan
            rates
    ***********************************************************************/
    tempfile boutique
    load_analysis_data boutique_rates
    save `boutique'

    * units inconsistency
    replace r = 100 * r if whichsource == "nonqual"

    /***********************************************************************
        Fill in bank and loan rates from earlier, then make graph
    ***********************************************************************/

    append using `bank'

    sort year whichsource grp

    keep if year == 2016

    gen barorder = cond(whichsource == "nonqual", _N - 1, ///
                   cond(whichsource == "savbnd", _N + 1, ///
                   cond(whichsource == "loan", _N, ///
                   cond(whichsource == "bank", _n, 0.5 + _n))))

    colorpalette "$f3", intensity(0.05(.05)1)
    local f31 = "`r(p8)'" 
    local f32 = "`r(p12)'"
    local f33 = "`r(p16)'"
    local f34 = "`r(p20)'"

    colorpalette "$u1", intensity(0.05(.05)1)
    local u11 = "`r(p8)'" 
    local u12 = "`r(p11)'"
    local u13 = "`r(p14)'"
    local u14 = "`r(p17)'"
    local u15 = "`r(p20)'"

    format r %9.2f

    #delimit ;
    twoway 
        (bar r barorder if whichsource == "bank" & grp == "0", 
            color("`f31'") barw(0.75) lc("$u4"))
        (bar r barorder if whichsource == "bank" & grp == "90", 
            color("`f32'") barw(0.75) lc("$u4"))
        (bar r barorder if whichsource == "bank" & grp == "99", 
            color("`f33'") barw(0.75) lc("$u4"))
        (bar r barorder if whichsource == "bank" & grp == "99.9", 
            color("`f34'") barw(0.75) lc("$u4"))
        (bar r barorder if whichsource == "boutique" & grp == "0", 
            color("`u11'") barw(0.75))
        (bar r barorder if whichsource == "boutique" & grp == "90", 
            color("`u12'") barw(0.75))
        (bar r barorder if whichsource == "boutique" & grp == "99", 
            color("`u13'") barw(0.75))
        (bar r barorder if whichsource == "boutique" & grp == "99.9", 
            color("`u14'") barw(0.75))
        (bar r barorder if whichsource == "boutique" & grp == "99.99", 
            color("`u15'") barw(0.75))
        (bar r barorder if whichsource == "nonqual", 
            color("$f1") lc("$f4") barw(0.75))
        (bar r barorder if whichsource == "savbnd", 
            color("$f3") lc("$f4") barw(0.75))
        (bar r barorder if whichsource == "loan", 
            color(white) lc("$f4") barw(0.75))
        (scatter r barorder, ms(none) mlab(r) mlabpos(12) mlabcolor(black))
        , 
        ytitle("Interest Rates (%)")
        ylab(, format(%9.0f))
        xtitle("")
        xlabel(0.5(2.5)12, nolabels notick)
        $gpr 
        legend(order(
                     1 "Bank, P0-90 Non-Int Wealth"
                     5 "Boutique, Bot 90 AGI"
                     9 "Boutique, Top 0.01 AGI"
                     2 "Bank, P90-99 Non-Int Wealth"
                     6 "Boutique, P90-99 AGI"
                     10 "Non-Qual Divs"  
                     3 "Bank, P99-99.9 Non-Int Wealth"
                     7 "Boutique, P99-99.9 AGI"
                     12 "Business Loans"
                     4 "Bank, Top 0.1 Non-Int Wealth"
                     8 "Boutique, P99.9-99.99 AGI"
                     11 "Savings Bonds")
           size(small) col(3) symxsize(*0.6) region(lc(white) margin(tiny)))
        xsize(6.5);
    #delimit cr
    graph_export "asset_interest_rates_boutique_breakout2016"

end/*%>*/

capture program drop graph_hetero_fixed/*%<*/
program define graph_hetero_fixed

    /***************************************************************************
        Average Rates of Return in Capitalized Data
    ***************************************************************************/

    /***********************************************************************
        Prep Alexi Savov deposit rates in e-mail received by EZ on 
            10/15/2019 with subject line "Re: Aggregate Deposits Data." 
            Collapse by mean to yield annual deposit rates, then collapse 
            again to get average of these. Store result in global macro.
    ***********************************************************************/
    local start_year 2016

    load_analysis_data deposit_rates
    keep if inrange(year, `start_year', 2016)
    
    collapse (mean) r_deposit // Then get average over years
    sum r_deposit
    local r_deposit = `r(mean)'

    /***********************************************************************
        Get averages of 10-year Treasury and Aaa yields over 
            1998-2016; store in global macro
    ***********************************************************************/

    load_analysis_data fred_rates
    keep if inrange(year, `start_year', 2016)

    collapse (mean) r_ust10 r_aaa

    sum r_ust10
    local r_ust10 = `r(mean)'

    sum r_aaa
    local r_aaa = `r(mean)'

    /***********************************************************************
        Prep tax data ranked by non-interest wealth and preferred wealth
    ***********************************************************************/

    load_analysis_data szz_niw
    keep if inrange(year, `start_year', 2016)

    keep if inlist(w${niw_defn_late}_group, 1, 7, 8, 9, 10, 5)
    assert inlist(group, "All", "P0-90", "P90-99", "P99-99.9", "P99.9-99.99", "P99.99-100")

    gen r_niw = (interest_info_scaled / taxbond_info) * 100

    gen groupnum = cond(w${niw_defn_late}_group != 1, w${niw_defn_late}_group - 5, ///
                                                      w${niw_defn_late}_group)

    replace groupnum = 6 if w${niw_defn_late}_group == 5

    keep groupnum r_niw

    tempfile inside_niw
    save `inside_niw'

    load_analysis_data szz
    keep if inrange(year, `start_year', 2016)

    keep if inlist(w${preferred_defn_late}_group, 1, 7, 8, 9, 10, 5)
    assert inlist(group, "All", "P0-90", "P90-99", "P99-99.9", "P99.9-99.99", "P99.99-100")

    gen r_pref = (interest_info_scaled / taxbond_info) * 100

    gen groupnum = cond(w${preferred_defn_late}_group != 1, w${preferred_defn_late}_group - 5, ///
                                                      w${preferred_defn_late}_group)

    * TODO: remove when niw file is fixed
    replace groupnum = 6 if w${preferred_defn_late}_group == 5


    keep groupnum r_pref

    tempfile inside_pref
    save `inside_pref'

    load_analysis_data szz_cmd
    keep if inrange(year, `start_year', 2016)

    keep if inlist(w${cmd_defn_late}_group, 1, 7, 8, 9, 10, 5)
    assert inlist(group, "All", "P0-90", "P90-99", "P99-99.9", "P99.9-99.99", "P99.99-100")

    gen r_cmd = ((fiint + intest) / taxbond_cmd_3tier) * 100

    gen groupnum = cond(w${cmd_defn_late}_group != 1, w${cmd_defn_late}_group - 5, ///
                                                      w${cmd_defn_late}_group)

    replace groupnum = 6 if w${cmd_defn_late}_group == 5

    keep groupnum r_cmd

    tempfile inside_cmd
    save `inside_cmd'

    /***********************************************************************
        Prep tax data ranked by AGI
    ***********************************************************************/
    load_analysis_data szz_agi
    keep if inrange(year, `start_year', 2016)

    keep year group interest_info_scaled taxbond_info

    * Groups are mutually exclusive, collectively exhaustive
    levelsof group, local(groups) clean
    assert "`groups'" == "P0-90 P90-99 P99-99.9 P99.9-99.99 P99.99-100"    

    /*******************************************************************
        Calculate total fixed income and fixed claims
    *******************************************************************/

    sort year
    foreach sumqty of varlist interest_info_scaled taxbond_info {
        egen total_`sumqty' = total(`sumqty')
    }

    /*************************************k*****************************
        Put totals in a dedicated row instead of Keeping
            them in a column
    *******************************************************************/

    expand 2 if group == "P0-90" 
    
    sort year group
    by year: assert group[1] == "P0-90" & group[2] == "P0-90"
    by year: replace group = "All" if _n == 1
    by year: replace interest_info_scaled = total_interest_info_scaled if group == "All"
    by year: replace taxbond_info = total_taxbond_info if group == "All"

    encode group, gen(groupnum)

    isid groupnum

    gen r_agi = (interest_info_scaled / taxbond_info) * 100

    keep groupnum r_agi

    tempfile inside_agi
    save `inside_agi'

    /***********************************************************************
        Merge deposit and SCF data files together, plus merge in 
            estate tax returns interest rates; then make plots
    ***********************************************************************/

    use `inside_pref', clear
    merge 1:1 groupnum using `inside_niw', assert(3) nogen
    merge 1:1 groupnum using `inside_cmd', assert(3) nogen
    merge 1:1 groupnum using `inside_agi', assert(3) nogen

    #delimit ;
    graph bar (asis) r_pref r_agi r_niw 
        ,
        $gpr
        bargap(20)
        bar(2, color("$p21"))
        bar(3, fcolor("$u8") lcolor(black))
        bar(1, color("$u7"))
        legend(order(1 "Baseline Wealth" 2 "AGI" 3 "Non-Interest Wealth" ) 
            region(lcolor(white) margin(tiny)) row(1)) 
        over(group, relabel(1 "All" 2 "P0-90" 3 "P90-99" 4 "P99-99.9" 5 "P99.9-99.99" 
                            6 "P99.99-100")) 
        ylab(0(1)5) yscale(range(0 5.3))
        ytitle("Return on Taxable-Interest Assets (%)")
        xsize(6.5);
    #delimit cr
    graph_export "avgfixreturn2016"

end/*%>*/

capture program drop graph_rates_overtime_fixed/*%<*/
program define graph_rates_overtime_fixed

    /***************************************************************************
        Interest rates
    ***************************************************************************/

    /***********************************************************************
        Prepare interest rates by wealth under our preferred 
            capitalization
    ***********************************************************************/

    load_analysis_data szz
    keep if inlist(w_group, 1, 7, 8, 3, 4, 5) 
    qui levelsof group, local(groups) clean
    assert "`groups'" == "All P0-90 P90-99 P99-100 P99.9-100 P99.99-100"

    #delimit ;
    replace group = cond(group == "All", "_all",
                    cond(group == "P0-90", "_bot90", 
                    cond(group == "P90-99", "_p90_99",
                    cond(group == "P99-100", "_top1", 
                    cond(group == "P99.9-100", "_top01", "_top001")))));
    #delimit cr

    gen taxbond_pref = cond(year < 2001, taxbond_cmd_3tier, taxbond_info)
    gen fixinc_pref = cond(year < 2001, inc_fix, interest_info_scaled)
    gen taxbond_eqmuf = taxbond + taxbond_mufmisc_sz

    keep year group fixinc_pref inc_fix taxbond_ini taxbond_pref taxbond_equal taxbond_eqmuf

    reshape wide inc_fix fixinc_pref taxbond_ini taxbond_pref taxbond_equal taxbond_eqmuf, i(year) j(group, string)

    assert inrange(fixinc_pref_all / inc_fix_all, 0.99, 1.01)
    drop inc_fix*

    gen fixinc_pref_bot99 = fixinc_pref_all - fixinc_pref_top1
    gen taxbond_pref_bot99 = taxbond_pref_all - taxbond_pref_top1

    foreach grp in top001 top01 top1 bot99 {
        gen r_`grp'_preferred = (fixinc_pref_`grp' / taxbond_pref_`grp') * 100
    }

    gen r_equreturns = (fixinc_pref_all / taxbond_pref_all) * 100
    gen r_equreturns_oldinttaxw = (fixinc_pref_all / taxbond_ini_all) * 100
    gen r_equreturns_sz20 = (fixinc_pref_all / taxbond_equal_all) * 100
    gen r_equreturns_mufmisc = (fixinc_pref_all / taxbond_eqmuf_all) * 100

    tempfile capitalized
    save `capitalized'

    /***********************************************************************
        Prep interest rates from FRED and Savov.
    ***********************************************************************/

    load_analysis_data fred_rates
    keep year r_aaa r_baa r_ust10 r_bbb
    keep if inrange(year, 1966, 2019)
    tempfile bondyields
    save `bondyields'

    load_analysis_data deposit_rates
    assert missing(r_deposit) if year == 1986
    drop if year == 1986
    tempfile deposit
    save `deposit'

    /***********************************************************************
        Merge data together and make graph(s)
    ***********************************************************************/

    use `bondyields', clear
    
    merge 1:1 year using `capitalized', assert(1 3) keepusing(r_*)
    assert year > 2016 if _merge == 1
    assert _merge == 3 if year < 2017
    drop _merge

    merge 1:1 year using `deposit', assert(1 3)
    assert !inrange(year, 1987, 2016) if _merge == 1
    assert _merge == 3 if inrange(year, 1987, 2016)
    drop _merge

    outsheet using interest_rate_series.csv, replace

    /***********************************************************************
        Make graph
    ***********************************************************************/

    sort year
    #delimit ;
    graph twoway 
        (connect r_equreturns_mufmisc year, color("$u3") ms(d))
        (connect r_ust10 year, ms(t) color("$u5") lp("-")) 
        (connect r_aaa year, ms(s) color("$u5") lp("-") lwidth(medthick)) 
        (connect r_baa year, c(l) ms(o) color("$u5") lp("-") lw(medthin)) 
        (connect r_deposit year, ms(x) mc("$f3") lc("$f3")) 
        (connect r_top001_preferred year, ms(o) color("$u1"))
        (connect r_top01_preferred year, ms(s) color("$u1"))
        (connect r_top1_preferred year, ms(oh) color("$u1") lp(_))
        (connect r_bot99_preferred year, ms(oh) color("$u1") lp("."))
            if year <= 2016,
        ytitle("Interest Rate (%)")
        xtitle(" ") xlab(1965(5)2015) ysca(range(-0.5 15.3)) ylabel(0(3)15) $gpr
        legend(order(6 "Top 0.01 Baseline" 7 "Top 0.1 Baseline" 
                    8 "Top 1 Baseline" 9 "P0-99 Baseline" 
                    1 "Equal Returns" 5 "Deposits" 
                    4 "Moody's Baa" 3 "Moody's Aaa" 2 "10-Yr. Treasury" 
                    )
            region(lcolor(white) margin(tiny)) row(3))
        xsize(6.5);
    #delimit cr
    graph_export interest_rate_series

end/*%>*/

capture program drop graph_cmd/*%<*/
program define graph_cmd

    /***************************************************************************
        CMD top and bottom rates with confidence interval
    ***************************************************************************/
    use $inputs/cmd_flows_rates_niwranks20210715.dta, clear

    keep year r_top01 r_bot99pt9 ub_r_top01 lb_r_top01 ub_r_bot99pt9 lb_r_bot99pt9
    rename (year r_top01 r_bot99pt9 ub_r_top01 lb_r_top01 ub_r_bot99pt9 lb_r_bot99pt9) ///
         (year r1 r2 ub_r1 lb_r1 ub_r2 lb_r2)

    sort year
    #delimit ;
    graph twoway (connect r1 year, msym(oh) mcolor("$u3") lcolor("$u3") lp(-))
        (connect ub_r1 year, msym(none) mcolor("$u3") lcolor("$u3") lp("."))
        (connect lb_r1 year, msym(none) mcolor("$u3") lcolor("$u3") lp("."))
        (connect r2 year, msym(sh) mcolor("$u1") lcolor("$u1") lp(-))
        (connect ub_r2 year, msym(none) mcolor("$u1") lcolor("$u1") lp("."))
        (connect lb_r2 year, msym(none) mcolor("$u1") lcolor("$u1") lp(".")), 
        $gpr
        ytitle("Interest Rate (%)") xtitle("") 
        xlab(1965(5)2015) xscale(range(1964 2017))
        yscale(range(-1 15)) ylab(0(3)15)
        legend(order(1 "Top 0.1% Non-Interest Wealth"
                     4 "Bottom 99.9% Non-Interest Wealth")
               region(lcolor(white) margin(tiny)) row(1))
        xsize(6.5);
    #delimit cr
    graph_export "cmd_interest_rates_with_bounds"

end/*%>*/

capture program drop graph_factors_overtime_fixed/*%<*/
program define graph_factors_overtime_fixed

    /***************************************************************************
        Capitalization factors
    ***************************************************************************/
    
    /***********************************************************************
        Calculate cap factors associated with CMD interest rates
    ***********************************************************************/

    use $inputs/cmd_flows_rates_niwranks20210715.dta, clear

    keep year r_top01 r_bot99pt9 ub_r_top01 lb_r_top01 ub_r_bot99pt9 lb_r_bot99pt9
    rename (year r_top01 r_bot99pt9 ub_r_top01 lb_r_top01 ub_r_bot99pt9 lb_r_bot99pt9) ///
         (year r1 r2 ub_r1 lb_r1 ub_r2 lb_r2)

    gen capfactor_cmd_top01niw = 100 / r1

    tempfile cmd
    save `cmd'

    /***********************************************************************
        Calculate equal returns, UST 10-year, and Moody's Aaa cap
            factors
    ***********************************************************************/

    import delimited using interest_rate_series.csv, clear // From 5C

    gen capfactor_equreturns = 100 / r_equreturns
    gen capfactor_equreturns_oldinttaxw = 100 / r_equreturns_oldinttaxw
    gen capfactor_equreturns_sz20 = 100 / r_equreturns_sz20
    gen capfactor_equreturns_mufmisc = 100 / r_equreturns_mufmisc
    gen capfactor_ust10 = 100 / r_ust10
    gen capfactor_aaa = 100 / r_aaa
    gen capfactor_baa = 100 / r_baa

    gen capfactor_top001_preferred = 100 / r_top001_preferred
    gen capfactor_top01_preferred = 100 / r_top01_preferred
    gen capfactor_top1_preferred = 100 / r_top1_preferred
    gen capfactor_bot99_preferred = 100 / r_bot99_preferred

    /***********************************************************************
        Merge CMD cap factors onto cap factors from above and make 
            graph
    ***********************************************************************/

    merge 1:1 year using `cmd', keepusing(capfactor_*) assert(1 3)
    assert year > 2016 if _merge == 1
    assert _merge == 3 if year <= 2016
    drop _merge

    outsheet year capfactor_* using fix_capfactors.csv, replace

    sort year
    #delimit ;
    twoway 
        (connect capfactor_equreturns_mufmisc year,  lc("$u3") mc("$u3") ms(d)) 
        (connect capfactor_ust10 year, ms(t) mc("$u5") lc("$u5") lp("-")) 
        (connect capfactor_aaa year, ms(s) lc("$u5") mc("$u5") lp("-") lw(medthick)) 
        (connect capfactor_baa year, ms(o) lc("$u5") mc("$u5") lp("-") lw(medthin)) 
        (connect capfactor_cmd_top01niw year, ms(oh) lc("$u1") mc("$u1") lpattern(".")) 
        (connect capfactor_top001_preferred year, ms(o) color("$u1"))
        (connect capfactor_top01_preferred year, ms(s) color("$u1") lp(-))
        (connect capfactor_top1_preferred year, ms(oh) color("$u1") lpattern(_))
            if inrange(year, 1966, 2016),
        $gpr
        ytitle("Capitalization Factor") ylab(0(25)125)
        xtitle("") xlab(1965(10)2015)
        legend(order(1 "Equal Returns"
                     8 "Top 1 Base"
                     7 "Top 0.1 Base" 
                     6 "Top 0.01 Base" 
                     5 "CMD Top 0.1"
                     4 "Moody's Baa" 3 "Moody's Aaa" 2 "10-Yr Treas" 
                    )
            region(lcolor(white) margin(tiny)) symxsize(*.5) col(4))
        xsize(6.5);
    #delimit cr
    graph_export "fix_capfactors"

end/*%>*/

capture program drop graph_return_ratio_fixed/*%<*/
program define graph_return_ratio_fixed

    /***************************************************************************
        Ratio of top group to macro returns
    ***************************************************************************/

    /***********************************************************************
        Retrieve macro rate and UST 10, Moody's Aaa from figure 5C 
            output CSV
    ***********************************************************************/

    import delimited using interest_rate_series.csv, clear

    rename r_equreturns r_macro
    drop if missing(r_macro)

    tempfile macrorate
    save `macrorate'

    /***********************************************************************
        Prepare CMD output
    ***********************************************************************/
    use $inputs/cmd_flows_rates_niwranks20210715.dta, clear

    keep year r_top01 r_bot99pt9 ub_r_top01 lb_r_top01 ub_r_bot99pt9 lb_r_bot99pt9
    rename (year r_top01 r_bot99pt9 ub_r_top01 lb_r_top01 ub_r_bot99pt9 lb_r_bot99pt9) ///
         (year r1 r2 ub_r1 lb_r1 ub_r2 lb_r2)
    isid year

    /***********************************************************************
        Merge ratios from different data sources and make plot
    ***********************************************************************/
    merge 1:1 year using `macrorate', keepusing(r_macro r_aaa r_baa r_ust10 r_top0*_preferred) ///
        assert(3) nogen

    foreach numerator of varlist r_aaa r_baa r_ust10 r1 ub_r1 lb_r1 r_top*_preferred {
        gen `numerator'_rmacro = `numerator' / r_macro
        assert `numerator'_rmacro > 0

        if regexm("`numerator'", "preferred") {
            replace `numerator'_rmacro = . if year < 2001 // Only want info returns

        }
    }

    gen revisionists_ratio = cond(year < 2003, 1, cond(year < 2008, 1.15, 1.4))

    #delimit ;
    twoway 
        (connect r_ust10_rmacro year, msym(t) lcolor("$u5") mcolor("$u5") lp("-"))
        (connect r_aaa_rmacro year, msym(s) lcolor("$u5") mcolor("$u5") lp("-") lw(medthick))
        (connect r_baa_rmacro year, msym(o) lcolor("$u5") mcolor("$u5") lp("-") lw(medthin))
        (connect r_top001_preferred_rmacro year, msym(o) mcolor("$u1") lcolor("$u1"))
        (connect r_top01_preferred_rmacro year, msym(s) mcolor("$u1") lcolor("$u1") lp("_"))
        (connect r1_rmacro year, msym(oh) mcolor("$u1") lcolor("$u1") lp("_"))
        (connect ub_r1_rmacro year, msym(none) lcolor("$u1") lpattern("."))
        (connect lb_r1_rmacro year, msym(none) lcolor("$u1") lpattern("."))
        ,
        yline(1, lcolor("$u3"))
        text(0.85 1973 "Equal returns assumption", size(small) color("$u3"))
        $gpr ytitle("Ratio of Rate to Equal Returns Rate") xtitle("") 
        xlab(1965(5)2020) xscale(range(1965 2020))
        yscale(range(-0.1 4)) ylabel(0(1)4)
        legend(order(3 "Moody's Baa" 4 "Top 0.01% Baseline" 
                     2 "Moody's Aaa" 5 "Top 0.1% Baseline"  
                     1 "10-year Treas."  6 "CMD Top 0.1% Non-Int Wealth") 
               region(lcolor(white)) row(3));
    #delimit cr
    graph_export "fix_top_to_macro_ratio"

end/*%>*/

capture program drop graph_comparison_trends/*%<*/
program define graph_comparison_trends

    tempfile imputation

    /***************************************************************************
        Effect of Refinement
    ***************************************************************************/

    /***********************************************************************
        Pull in collapsed tax data and rename/generate SZ and SZZ
            measures for housing wealth, pass-through wealth, and pension
            wealth. Drop variables which are not important for comparison.
    ***********************************************************************/
    load_analysis_data szz
    keep if inlist(group, "All", "P99.9-100")
    replace group = cond(group == "All", "_all", "_top01")

    gen taxbond_pref = cond(year < 2001, ///
                               taxbond_cmd_3tier, ///
                               taxbond_info)

    gen taxbond_eqmuf = taxbond + taxbond_mufmisc_sz

    keep year group hweal20 hweal_preferred taxbond taxbond_ini taxbond_ust10 taxbond_baa taxbond_aaa ///
        taxbond_pref taxbond_cmd_3tier taxbond_cmd_ub taxbond_cmd_lb ///
        taxbond_equal taxbond_eqmuf

    rename taxbond_cmd_3tier taxbond_cmd3

    reshape wide hweal20 hweal_preferred taxbond taxbond_ini taxbond_ust10 taxbond_baa taxbond_aaa ///
        taxbond_pref taxbond_cmd3 taxbond_cmd_ub taxbond_cmd_lb ///
        taxbond_equal taxbond_eqmuf, i(year) j(group, string)

    /***********************************************************************
       Compute shares of net household with attributable to fixed
            income and c-corp eqåuity wealth categories among top 0.1% under
            different capitalization assumptions
    ***********************************************************************/

    foreach taxbond in taxbond taxbond_ini taxbond_pref taxbond_baa taxbond_aaa taxbond_ust10 ///
        taxbond_cmd3 taxbond_cmd_ub taxbond_cmd_lb taxbond_equal taxbond_eqmuf {
        gen top01_`taxbond'_sh_aggwlth = `taxbond'_top01 / hweal_preferred_all * 100
    }

    replace top01_taxbond_ini_sh_aggwlth = taxbond_ini_top01 / hweal20_all * 100

    keep year top01_*_sh_aggwlth
    save `imputation'

    /************************************************************************
        Make graphs
    ************************************************************************/

    sort year
    #delimit ;
    twoway 
           (connect top01_taxbond_eqmuf_sh_aggwlth year, lc("$u3") mc("$u3") ms(s)) 
           (connect top01_taxbond_ust10_sh_aggwlth year, ms(t) lp("-") mc("$u5") lc("$u5")) 
           (connect top01_taxbond_aaa_sh_aggwlth year, ms(s) lp("-") mc("$u5") lc("$u5") lw(medthick)) 
           (connect top01_taxbond_baa_sh_aggwlth year, ms(o) lp("-") mc("$u5") lc("$u5") lw(medthin)) 
           (connect top01_taxbond_pref_sh_aggwlth year if year >= 2001, ms(O) lcolor("$u1") mcolor("$u1")
                lwidth(medthick))
           (connect top01_taxbond_cmd3_sh_aggwlth year, ms(oh) lcolor("$u1") mcolor("$u1") lp("."))
           (connect top01_taxbond_cmd_ub_sh_aggwlth year, ms(sh) lcolor("$u1") mcolor("$u1") lp("-"))
           (connect top01_taxbond_cmd_lb_sh_aggwlth year, ms(sh) lcolor("$u1") mcolor("$u1") lp("-"))
                ,
        ytitle("Top 0.1% Taxable Fixed Income Wealth as" "Share of Net Household Wealth (%)") 
        xtitle(" ") 
        xlab(1965(5)2015, labsize(small)) xsca(range(1964 2016))
        ylab(0(1)7) ysca(range(-0.5 8))
        legend(region(lcolor(white) margin(tiny)) row(3) symxsize(*.5)
            order(2 "10-Yr. Treas" 1 "Equal Returns" 7 "CMD 2-Tier Hi" 
                  3 "Moody's Aaa" 5 "Info Returns"  8 "CMD 2-Tier Low"
                  4 "Moody's Baa" 6 "CMD 3-Tier")) 
        $gpr;
    #delimit cr
    graph_export "effect_refinement_fix"

end/*%>*/

capture program drop graph_mse/*%<*/
program define graph_mse 
    
    /***************************************************************************
        Predictions and error in the SCF
    ***************************************************************************/

    /***********************************************************************
        Get total adults from parameters
    ***********************************************************************/

    import excel using $inputs/parameters.xlsx, ///
        firstrow clear

    drop if missing(yr)
    rename yr year
    keep year totadults20

    replace totadults20 = totadults20 * 1E3 // Scale out of thousands

    tempfile parameters
    save `parameters'

    /***********************************************************************
        Get top capitalization factors from figure 6.A output
    ***********************************************************************/

    import delimited using fix_capfactors.csv, clear

    keep if inrange(year, 1998, 2016)

    gen sz20capfactor_top1 = cond(year < 2003, capfactor_equreturns, ///
                              cond(year < 2008, capfactor_equreturns / 1.15, capfactor_equreturns / 1.4))

    keep year sz20capfactor_top1 capfactor_cmd_top01niw capfactor_equreturns ///
        capfactor_equreturns_oldinttaxw capfactor_equreturns_mufmisc

    tempfile maincapfactors
    save `maincapfactors'

    /***********************************************************************
        Calculate bottom CMD capitalization factor: this is a 
            residual after capitalizing the top 0.1% ranked by non-interest
            wealth
    ***********************************************************************/

    load_taxdata, rankspec($niw_defn_late) startyr(2001) endyr(2016)

    tempfile niw_late
    save `niw_late'

    load_taxdata, rankspec($niw_defn_mid) startyr(1998) endyr(2000)

    append using `niw_late'

    gen w_group = cond(year < 2001, w${niw_defn_mid}_group, w${niw_defn_late}_group)
    assert !missing(w_group)
    keep if inlist(w_group, 1, 4)
    assert inlist(group, "All", "P99.9-100")

    gen fixinc = fiint + intest

    keep year w_group taxbond_cmd fixinc

    reshape wide taxbond_cmd fixinc, i(year) j(w_group)
    rename (*1 *4) (*_all *_top01)

    gen capfactor_cmd_bot99pt9niw = (taxbond_cmd_all - taxbond_cmd_top01) / (fixinc_all - fixinc_top01)

    keep year capfactor_cmd_bot99pt9niw

    tempfile capfactor_cmd_bot99pt9niw
    save `capfactor_cmd_bot99pt9niw'

    /***********************************************************************
        Implement capitalization in SCF
    ***********************************************************************/

    /*******************************************************************
        Load data, restrict to years of interest, and get
            interest and wealth ranks
    *******************************************************************/

    use $inputs/scf_revision.dta, clear

    keep if inrange(year, 1998, 2016)

    es_rank_scf, rankvar(networth_pref) othersplitvars(intinc inttaxw) outname(es_wlthrank)

    cumul niw [aw = wgt], by(year) gen(niw_rank_baseline)

    merge m:1 year using `parameters', keepusing(totadults20) assert(2 3) keep(3) nogen

    sort year
    qui egen num_ind = total(wgt), by(year) // Total individuals according to SCF 

    qui gen pct_ind_above = 1 - niw_rank_baseline

    qui gen num_ind_above = pct_ind_above * num_ind

    qui gen pct_ind_above_taxdata = num_ind_above / (totadults20 - 400)

    qui gen es_niwrank = 1 - pct_ind_above_taxdata
    assert inrange(es_niwrank, 0, 1)

    /*******************************************************************
        Apply different capitalizations
    *******************************************************************/

    foreach mergefile in maincapfactors capfactor_cmd_bot99pt9niw {
        merge m:1 year using ``mergefile'', assert(2 3)
        assert mod(year - 1989, 3) > 0 if _merge == 2
        assert _merge == 3 if mod(year - 1989, 3) == 0 
        keep if _merge == 3
        drop _merge
    }

    * For this plot we only need top 1%ers and up in terms of wealth
    drop if es_wlthrank < 0.99

    gen taxbond_cmd = cond(es_niwrank >= 0.999, intinc * capfactor_cmd_top01niw, ///
                           intinc * capfactor_cmd_bot99pt9niw)

    gen taxbond_equret = intinc * capfactor_equreturns

    gen taxbond_psz18 = intinc * capfactor_equreturns_oldinttaxw

    gen taxbond_sz20 = intinc * sz20capfactor_top1

    gen taxbond_mufmisc = intinc * capfactor_equreturns_mufmisc

    tempfile orig
    save `orig'

    /*******************************************************************
        Collapse to yield aggregates by interest and wealth 
            groups
    *******************************************************************/

    foreach wlthgrp in top1 top01 top001 {

        local threshold = cond("`wlthgrp'" == "top1", 0.99, cond("`wlthgrp'" == "top01", 0.999, 0.9999))

        use `orig', clear

        if "`wlthgrp'" == "top1" {
            assert es_wlthrank >= `threshold'
        }
        else {
            drop if es_wlthrank < `threshold'
        }

        #delimit ;
        collapse (mean) taxbond_actual_`wlthgrp' = inttaxw
            taxbond_cmd_`wlthgrp' = taxbond_cmd 
            taxbond_psz18_`wlthgrp' = taxbond_psz18
            taxbond_equret_`wlthgrp' = taxbond_equret
            taxbond_sz20_`wlthgrp' = taxbond_sz20 
            taxbond_mufmisc_`wlthgrp' = taxbond_mufmisc [aw = wgt], by(year);
        #delimit cr

        foreach intomillions of varlist taxbond_* {
            replace `intomillions' = `intomillions' / 1E6
        }

        tempfile `wlthgrp'
        save ``wlthgrp''
    }

    use `top1', clear
    merge 1:1 year using `top01', assert(3) nogen
    merge 1:1 year using `top001', assert(3) nogen

    /***********************************************************************
        Make plots
    ***********************************************************************/

    #delimit ;
    twoway  (function y=x,  range(0 7) lcolor(gs10) lpattern(dash))
        (scatter taxbond_cmd_top1 taxbond_actual_top1, 
            ms(oh) color("$u1")) 
        (scatter taxbond_mufmisc_top1 taxbond_actual_top1, 
             ms(sh) color("$u3"))
        (scatter taxbond_sz20_top1 taxbond_actual_top1, 
            ms(x) color("$p2"))
        (scatter taxbond_cmd_top01 taxbond_actual_top01, 
            ms(d) color("$u1")) 
        (scatter taxbond_mufmisc_top01 taxbond_actual_top01, 
            ms(t) color("$u3"))
        (scatter taxbond_sz20_top01 taxbond_actual_top01, 
            ms(+) color("$p2"))
        (scatter taxbond_cmd_top001 taxbond_actual_top001, 
            ms(o) color("$u1") mlab(year) mlabcol("$u1") mlabpos(6))
        (scatter taxbond_mufmisc_top001 taxbond_actual_top001, 
            ms(s) color("$u3") mlab(year) mlabcol("$u3") mlabpos(12))
        (scatter taxbond_sz20_top001 taxbond_actual_top001, 
            ms(X) color("$p2") mlab(year) mlabcol("$p2") mlabpos(12))
        ,
        xtitle("Actual Per-Capita Mean (Millions)")
        ytitle("Predicted Per-Capita Mean (Millions)")
        graphregion(lcolor(white) fcolor(white)) plotregion(color(white))
        xlab(0(1)7) 
        legend(order(
                3 "P99-100 Equal Returns"   6 "P99.9-100 Equal Returns"  9 "P99.99-100 Equal Returns"
                4 "P99-100 SZ20"         7 "P99.9-100 SZ20"        10 "P99.99-100 SZ20"
                2 "P99-100 CMD 2-Tier"      5 "P99.9-100 CMD 2-Tier"      8 "P99.99-100 CMD 2-Tier"
                ) 
            region(lcolor(white)) row(4) size(small));
    #delimit cr
    graph_export fix_predictions_wealthranks

end/*%>*/

/*%>*/
*******************************************************************************

****************************************************************************
* Pass-through
****************************************************************************
capture program drop graph_aggregate_passthru/*%<*/
program define graph_aggregate_passthru

    /***************************************************************************
        Aggregate pass-through business across data sources
    ***************************************************************************/
    
    /***********************************************************************
        Pull PSZ 2018 parameter aggregates for S-corporations, 
            C-corporations, and non-corporate business. Assume 20% of
            C-corporations are private. 
    ***********************************************************************/
    load_analysis_data parameters_older
    keep year ttdivw ttscorw ttschcpartw
    keep if inrange(year, 1989, 2015)
    rename (ttdivw ttscorw ttschcpartw) (ttdivw_psz18 ttscorw_psz18 ttschcpartw_psz18)

    tempfile psz18parameters
    save `psz18parameters'

    load_analysis_data parameters_old

    keep if inrange(year, 1989, 2019)

    keep year ttscorw ttschcpartw 

    tempfile parameters
    save `parameters'

    /***********************************************************************
        Pull USFA market value of closely held C-corporations; 
            ttdivw encompassses all C-corporation equity, which is too 
            broad
    ***********************************************************************/
    load_analysis_data usfa

    keep if year >= 1989

    tempfile usfa
    save `usfa'

    /***********************************************************************
        Pull SCF aggregates
    ***********************************************************************/
    load_analysis_data scfrevision

    collapse (sum) hwbus_scf = hwbus pthru_scf = pthru privccorw_scf = privccorw ///
        [fw = wgt1B], by(year)

    foreach billionX of varlist *_scf { // Scale down billion-weighted aggregates
        replace `billionX' = `billionX' / 1E9
    }

    * Affirm that only diff b/w pthru and hwbus is inclusion of private C-corps
    gen privccorw_check = hwbus_scf - pthru_scf
    assert inrange(privccorw_check / privccorw_scf, 0.999, 1.001)
    drop privccorw_*

    tempfile scf 
    save `scf'

    /***********************************************************************
        Get totals with SZZ data
    ***********************************************************************/
    load_analysis_data szz 

    keep if group == "All" & year >= 1989

    assert missing(s_value_avg_ebitda_red) & missing(p_value_avg_ebitda_red) if year < 2001

    gen pthru_szz = s_value_avg_ebitda_hybnof_red + p_value_avg_ebitda_hybnof_red + ///
        solepropw + missing_scorp + missing_pship 

    gen pthru_szz_check = hwbus_pref - 0.2 * ccorw_pref
    assert inrange(pthru_szz / pthru_szz_check, 0.99, 1.01) if pthru_szz != .
    replace pthru_szz = hwbus_pref - 0.2 * ccorw_pref if year < 2001

    rename hwbus_pref hwbus_szz

    * Make a version w/o liquidity discount
    gen pthru_szz_nodisct = (1 / 0.9) * (s_value_avg_ebitda_full_red + p_value_avg_ebitda_full_red) + ///
                            solepropw + (1 / 0.9) * (missing_scorp + missing_pship)

    gen pthru_szz_half = (s_value_avg_ebitda_half_red + p_value_avg_ebitda_half_red) + ///
                            solepropw + (missing_scorp + missing_pship)

    gen hwbus_szz_nodisct = pthru_szz_nodisct + (0.2 * ccorw_pref)

    gen pthru_szz_hyb = s_value_avg_ebitda_hyb_red + p_value_avg_ebitda_hyb_red + ///
        solepropw + missing_scorp + missing_pship 

    keep year hwbus_szz* pthru_szz* hwbus_base* pthru_base*

    tempfile szz
    save `szz'

    /***********************************************************************
        Load aggregates scaled by national income from figure 3 to 
            affirm that this matches up
    ***********************************************************************/
    import delimited aggwealth.csv, clear

    tempfile figure3series
    save `figure3series'

    /***********************************************************************
        Merge data together
    ***********************************************************************/
    load_analysis_data natinc
    tempfile natinc 
    replace y = 1e-12 * y
    save `natinc'

    use `parameters', clear
    
    merge 1:1 year using `psz18parameters', assert(1 3)
    assert year > 2015 if _merge == 1
    assert _merge == 3 if year <= 2015
    drop _merge
    
    merge 1:1 year using `usfa', assert(3) nogen

    gen pthru_psz18 = ttscorw_psz18 + ttschcpartw_psz18
    gen hwbus_psz18 = pthru_psz18 + privccorw_usfa 
    drop tt*_psz18

    gen pthru_parameters = ttscorw + ttschcpartw
    gen hwbus_parameters = pthru_parameters + privccorw_usfa
    drop ttscorw ttschcpartw privccorw_usfa

    merge 1:1 year using `szz', assert(1 3) nogen

    merge 1:1 year using `scf', assert(1 3)
    assert mod(year - 1989, 3) > 0 if _merge == 1
    assert _merge == 3 if mod(year - 1989, 3) == 0
    drop _merge

    merge 1:1 year using `natinc', keep(3) nogen

    /***********************************************************************
        Make graph
    ***********************************************************************/

    foreach millions of varlist *_parameters *_psz18 {
        replace `millions' = `millions' / 1E6 // currently in M, scale to T
    }

    * currently in dollars, scale to T
    foreach dollars of varlist *_scf *_szz *_szz_nodisct *_szz_half *_szz_hyb *_base {
        replace `dollars' = `dollars' / 1E12
    }

    * As a share of national income
    foreach scalebyni of varlist pthru_* hwbus_* {
        gen `scalebyni'_sh_ni = (`scalebyni' / y) * 100
    }        

    /***********************************************************************
        Ensure national income share totals match what we showed in
            figure 3
    ***********************************************************************/

    merge 1:1 year using `figure3series', keep(1 3) keepusing(ttbus_share_y ttbus_preferred_share_y) nogen 

    gen preferred_check = pthru_base_sh_ni / ttbus_preferred_share_y
    gen sz2020_check = pthru_parameters_sh_ni / ttbus_share_y

    assert inrange(preferred_check, 0.99, 1.01) & inrange(sz2020_check, 0.99, 1.01) if inrange(year, 2001, 2016)
    drop ttbus_preferred_share_y ttbus_share_y *_check

    gen year_m15 = year - .15
    gen year_p15 = year + .15

    gen year_m3 = year - 0.3
    gen year_p3 = year + 0.3

    sort year
    #delimit ;
    graph twoway (bar pthru_scf_sh_ni year_m3 if mod(year - 1989, 3) == 0 & year < 2016, 
            color("$p21") barw(.3))
        (bar pthru_parameters_sh_ni year if mod(year - 1989, 3) == 0 & year < 2016, 
            fcolor("$u8") lcolor(black) barw(.3))
        (bar pthru_parameters_sh_ni year_m15 if mod(year - 1989, 3) != 0 & year < 2016, 
            fcolor("$u8") lcolor(black) barw(.3))
        (bar pthru_scf_sh_ni year_m15 if mod(year - 1989, 3) == 0 & year >= 2016, 
            color("$p21") barw(.3))
        (bar pthru_parameters_sh_ni year_p15 if year >= 2016, 
            fcolor("$u8") lcolor(black) barw(.3))
        (connect pthru_szz_sh_ni year if year >= 2000, lc("$u1") lw(thick) mc("$u1") ms(O))
        (connect pthru_szz_sh_ni year if year < 2001, lc("$u1") lw(thick) mc("$u1") lp(-) ms(S))
        (connect pthru_szz_nodisct_sh_ni year, color("$u1") ms(Oh) lp(_))
        (connect pthru_szz_hyb_sh_ni year, color("$u1") ms(Dh) lp(.-)) 
        ,
        $gpr 
        xtitle("")
        ytitle("Aggregate Value / National Income (%)") 
        xsca(range(1988 2020))
        xlab(1989(3)2019)
        ylab(0(20)120)
        legend(order(1 "Pass-through and Informal Business (SCF)" 
                     2 "Proprietorships, Partnerships, S-corps (USFA/Baseline/SZ20)" 
                     6 "Proprietorships; Supp. Pships + S-corps, Info Returns" 
                     7 "Proprietorships; Supp. Pships + S-corps w/o Info Returns" 
                     9 "Proprietorships; Supp. Pships + S-corps w/o Fin Adj"
                     8 "Proprietorships; Supp. Pships + S-corps w/o Liquidity, Fin, Labor Adj"
                     )
                region(lcolor(white) margin(tiny)) row(3) symxsize(*0.5)) 
        xsize(10);  
    #delimit cr
    graph_export "privbiz_aggregates"

end/*%>*/

capture program drop graph_indy_passthru/*%<*/
program define graph_indy_passthru
	xx
    load_analysis_data hetero_vals_v3

    local liq = .1

    * To harmonize units with SOI collapse
    replace profits_s = 1e-3 * profits_s

    * Compute return on equity by industry
    gen return_on_equity_avg = 1e-6 * profits_S_soi/( (1-`liq') * value_s_modelaverage)
    impose_bounds return_on_equity_avg, upper(1) lower(0)

    gen return_on_equity = 1e-6 * profits_S_soi/( (1-`liq') * value_s_ebitd)
    impose_bounds return_on_equity, upper(1) lower(0)

    * Compute aggregate return on equity
    bys year: egen total_profits = total(profits_S_soi)

    bys year: egen total_value_s_avg = total(value_s_modelaverage)
    gen return_on_equity_agg_avg = 1e-6 * total_profits/( (1-`liq') * total_value_s_avg)

    bys year: egen total_value_s = total(value_s_ebitd)
    gen return_on_equity_agg = 1e-6 * total_profits/( (1-`liq') * total_value_s)

    tempfile scorps
    save `scorps'

    use `scorps', clear
    keep if year >= 2001 & year <= 2016
    usd2016, var(value_s_modelaverage) yr(year)
    collapse (mean) value_s_modelaverage return_on_equity_avg return_on_equity_agg_avg, ///
        by(naics_4d longname)
    * Dropping mgmt and hold co
    drop if naics_4d == 5511

    * graph: compare to hetero return, perhaps using thin bar graph with only
    * top 40 industries by aggregate market cap

    replace return_on_equity_avg = 100 * return_on_equity_avg
    replace return_on_equity_agg_avg = 100 * return_on_equity_agg_avg
    gsort - value_s_modelaverage
    sum return_on_equity_agg_avg
    local roe = `r(mean)'

    #delimit ;
    graph hbar return_on_equity_avg if _n <= 25, 
        over(longname, sort(return_on_equity_avg) reverse label(labs(medsmall))) 
        $gpr xsize(2.8) ylab(0(5)35)
        yline(`roe', lc("$u3") lw(thick))
        bar(1, color("$u1"))
        ytitle("S-corp Returns (%)");
    #delimit cr
	graph_export "hetero_return_scorp_indy"
    

end/*%>*/

capture program drop graph_passthru_hockey/*%<*/
program define graph_passthru_hockey
    
    load_analysis_data pthru_returns

    bys groupvar: egen total_pthru_val = total(pthru_val)
    gen pthru_val_share = 100 * pthru_val / total_pthru_val

    gsort groupvar group

    replace pthru_returns = 100 * pthru_returns

    #delimit ;
    twoway
        (connect pthru_returns group if groupvar == "wpthru_inc", 
            ms(small) mc("$u3") lc("$u3") lw(thin) ms(s))
        (connect pthru_returns group if groupvar == "wagi",
            ms(small) mc("$p2") lc("$p2") lw(thin) ms(d)) 
        (connect pthru_returns group if groupvar == "whweal112", 
            ms(small) mc("$u1") lc("$u1") lw(thin) ms(o)) 
        if group > 75, 
        xtitle("Percentile") ytitle("Returns (Pass-Through Income / Value, %)")
        ysca(range(-1 25)) ylab(0(3)24)
        xsca(range(74 101))
        legend(order(3 "Baseline Wealth" 2 "AGI" 1 "Pass-Through Income") col(3)
        region(lc(white))) $gpr xsize(5.5);
    graph_export "passthru_returns_ranked"

    * Minor issue with imprecision in these groups => Be careful with group2 def;
    gen group2 = cond(group > 99.95, 5,
                 cond(group > 99.85, 4,  
                 cond(group > 98, 3, 
                 cond(group > 89, 2, 
                 cond(group > 0, 1, 0)))));

    collapse (sum) pthru_val_share, by(group2 groupvar);

    reshape wide pthru_val_share, i(group2) j(groupvar, string);

    graph bar (asis) 
        pthru_val_sharewhweal112 
        pthru_val_sharewagi
        pthru_val_sharewpthru_inc
        ,
        $gpr
        bargap(20)
        bar(1, color("$u7")) 
        bar(2, color("$p21")) 
        bar(3, fcolor("$u8") lcolor(black)) 
        legend(order(1 "Baseline Wealth" 2 "AGI" 3 "Pass-Through Income") 
            region(lcolor(white) margin(tiny)) row(1)) 
        over(group2, relabel(1 "P0" 2 "P1-90" 3 "P90-99" 
                            4 "P99-99.9" 5 "P99.9-99.99" 6 "P99.99-100")) 
        ylab(0(5)30) yscale(range(0 31)) 
        ytitle("Share of Pass-Through Value (%)")
        xsize(5.5);
    #delimit cr

    graph_export "passthru_wealthshare_ranked"
	
end/*%>*/


*******************************************************************************
* Results and Comparisons%<
*******************************************************************************
capture program drop table_topwealth/*%<*/
program define table_topwealth

    /*******************************************************************************
        Thresholds and Average Wealth in Top Wealth Groups

            Cycles through wealth specifications (Saez-Zucman baseline, Saez-Zucman
            with our private business adjustments, and our preferred measures) in 
            2016 to make table 1 according to a variety of specifications. In 
            original draft, table reports 2014 values.
    *******************************************************************************/

    if "`1'" == "" {
        local year 2016
    }
    else {
        local year `1'
    }

    /***************************************************************************
        Prep to build head of the table (e.g. title, etc.) via listtex, 
            which doesn't depend on data
    ***************************************************************************/

    clear
    set obs 1

    gen tab = "\begin{tabular}{llllllcccc}"
    gen hline = "\hline"
    gen toprule = "\toprule"
    gen midrule = "\midrule"
    gen botrule = "\bottomrule"
    gen cline = "\cline{4-6} \cline{8-10}"
    gen cmidrule = "\cmidrule(lr){4-6} \cmidrule(lr){8-10}"

    gen c = " "

    #delimit ;
    gen title1 = "Wealth group & Count & Baseline Threshold & \multicolumn{3}{c}{Average wealth} 
                    & &  \multicolumn{3}{c}{Wealth share}";
    #delimit cr
    gen title2 = "& & & Baseline & Equal Returns & Revised SZ & & Baseline & Equal Returns & Revised SZ"

    local tabname = "thresholdsavgwealth`year'"

    *Output
    listtex tab if _n == 1 using "`tabname'.tex", replace rstyle(none)
    listtex toprule if _n == 1, appendto("`tabname'.tex") rstyle(none)
    listtex title1 if _n == 1, appendto("`tabname'.tex") rstyle(tabular)
    listtex cmidrule if _n == 1, appendto("`tabname'.tex") rstyle(none)
    listtex title2 if _n == 1, appendto("`tabname'.tex") rstyle(tabular)
    listtex midrule if _n==1 , appendto("`tabname'.tex") rstyle(none)    

    /***************************************************************************
        Cycling through ranking specifications, make data-dependent 
            parts of table, namely wealth totals, levels, thresholds, and group
            counts under each ranking specification: SZ (baseline) and preferred
    ***************************************************************************/

    tempfile hweal$preferred_defn_late hweal$sz_wlth_defn_late_v3
    foreach rankspec in $preferred_defn_late $sz_wlth_defn_late_v3 {

        /***********************************************************************
            Load in tax data using dina_collapse program for year 2016.
                Ensure 10 wealth groups. 
        ***********************************************************************/

        load_taxdata, rankspec(`rankspec') year(`year')

        assert _N == 11

        /***********************************************************************
            Drop unnecessary variables
        ***********************************************************************/

        keep group w`rankspec'_group n threshold hweal$sz_wlth_defn_late_v3 ///
            hweal$preferred_defn_late

        /***********************************************************************
            Compute wealth shares and average wealth in groups
        ***********************************************************************/

        sum n if w`rankspec'_group == 1
        local aggn = `r(mean)'
        sum hweal$preferred_defn_late if w`rankspec'_group == 1
        local aggwlth = `r(mean)'

        * Note: after this step, hweal`spec' is A MEAN and NOT A TOTAL
        foreach hweal of varlist hweal* {
            gen sh_`hweal'  = 100 * (`hweal' / `hweal'[1])
            replace `hweal' = (`hweal' / n)
        }

        /***********************************************************************
            Reformat variables in preparation to put in table
        ***********************************************************************/

        replace threshold = cond(group == "All", ., round(threshold, 10^3))

        sort w`rankspec'_group

        forv i = 2 / 5 {
            local j = `i' + 6
            replace threshold = threshold[`i'] in `j'
        }

        * Label wealth groups
        gen lab = cond(group == "All", "Full population", ///
                    cond(group == "P90-100", "Top 10\%", ///
                    cond(group == "P99-100", "Top 1\%", ///
                    cond(group == "P99.9-100", "Top 0.1\%", ///
                    cond(group == "P99.99-100", "Top 0.01\%", ///
                    cond(group == "P99.999-100", "Top 0.001\%", ///
                    cond(group == "P0-90", "Bottom 90\%", ///
                    cond(group == "P90-99", "Top 10-1\%", ///
                    cond(group == "P99-99.9", "Top 1-0.1\%", ///
                    cond(group == "P99.9-99.99", "Top 0.1-0.01\%", ///
                    "Top 0.01-0.001\%"))))))))))
        
        foreach hweal of varlist hweal* {
            recast double `hweal'
            replace `hweal' = round(`hweal', 1000)
        }

        rename w`rankspec'_group w_group

        * formatting
        gen n_raw = n
        recast double n
        replace n = round(n, 100)

        format threshold* hweal* n* %15.0fc
        format sh* %15.1f

        /* Converting to strings, setting blanks properly, and adding dollar signs 
            in LaTeX-readable format */
        tostring threshold* hweal* sh*, replace usedisplayformat force
        replace threshold = "" if threshold == "."

        foreach dollarvar of varlist threshold hweal* {
            replace `dollarvar' = " \\$" + `dollarvar' if `dollarvar' != ""
        }

        replace sh_hweal$sz_wlth_defn_late_v3 = sh_hweal$sz_wlth_defn_late_v3 + "\%"
        replace sh_hweal$preferred_defn_late = sh_hweal$preferred_defn_late + "\%"

        save `hweal`rankspec'', replace

    }

    * Merge on SZ2020 data and compute extra columns 
    tempfile sz2020
    load_analysis_data sz2020_top
    keep if year == 2016
    expand 11

    gen w_group = _n
    rename *_pszes20 *

    * Top groups
    gen sh_sz2020 = 100 if _n == 1
    replace sh_sz2020 = top10 if _n == 2
    replace sh_sz2020 = top1 if _n == 3
    replace sh_sz2020 = top01 if _n == 4
    replace sh_sz2020 = top001 if _n == 5
    replace sh_sz2020 = top0001 if _n == 6

    * Intermediates
    replace sh_sz2020 = 100 - top10 if _n == 7
    replace sh_sz2020 = top10 - top1 if _n == 8
    replace sh_sz2020 = top1 - top01 if _n == 9
    replace sh_sz2020 = top01 - top001 if _n == 10
    replace sh_sz2020 = top001 - top0001 if _n == 11

    keep year w_group sh_sz2020
    save `sz2020'

    use `hweal$preferred_defn_late', clear
    drop *$sz_wlth_defn_late_v3
    merge 1:1 w_group using `hweal$sz_wlth_defn_late_v3', keep(3) ///
        keepusing(*$sz_wlth_defn_late_v3) nogen
    merge 1:1 w_group using `sz2020', keep(1 3) nogen

    gen hweal_sz2020 = 1e-2 * ( sh_sz2020 * `aggwlth' / n_raw)
    recast double hweal_sz2020
    replace hweal_sz2020 = round(hweal_sz2020, 1000)

    format hweal_sz2020 %15.0fc
    format sh_sz2020 %15.1f

    tostring hweal_sz2020 sh_sz2020, replace usedisplayformat force
    replace threshold = "" if threshold == "."

    foreach dollarvar of varlist hweal_sz2020 {
        replace `dollarvar' = " \\$" + `dollarvar' if `dollarvar' != ""
    }

    replace sh_sz2020 = sh_sz2020 + "\%"

    /***********************************************************************
        Write data-dependent part of table under ranking in question
            via listtex
    ***********************************************************************/

    gen panel1 = "\multicolumn{10}{l}{\textbf{Panel A. Top wealth groups}}" in 1
    gen panel2 = "\multicolumn{10}{l}{\textbf{Panel B. Intermediate wealth groups}}" in 1

    gen c = " " in 1
    gen vsp = "\vspace{-.1in}" in 1

    * Panel 1 
    listtex panel1 if _n == 1 , appendto("`tabname'.tex") rstyle(none)
    listtex c if _n == 1, appendto("`tabname'.tex") rstyle(tabular) 
    listtex lab n threshold hweal$preferred_defn_late hweal$sz_wlth_defn_late_v3 ///
        hweal_sz2020 c sh_hweal$preferred_defn_late sh_hweal$sz_wlth_defn_late_v3 ///
        sh_sz2020 if w_group <= 6, ///
            appendto("`tabname'.tex") rstyle(tabular)

    * Panel 2
    listtex vsp if _n == 1, appendto("`tabname'.tex") rstyle(tabular)   
    listtex panel2 if _n == 1 , appendto("`tabname'.tex") rstyle(tabular)
    listtex c if _n == 1, appendto("`tabname'.tex") rstyle(none) 
    listtex lab n threshold hweal$preferred_defn_late hweal$sz_wlth_defn_late_v3 ///
        hweal_sz2020 c sh_hweal$preferred_defn_late sh_hweal$sz_wlth_defn_late_v3 ///
        sh_sz2020 if w_group > 6, ///
            appendto("`tabname'.tex") rstyle(tabular)

    /***************************************************************************
        Close out table
    ***************************************************************************/

    gen hline = "\hline" in 1
    gen botrule = "\bottomrule"
    gen end = "\end{tabular}" in 1

    listtex botrule if _n == 1, appendto("`tabname'.tex") rstyle(none)  
    listtex end if _n == 1, appendto("`tabname'.tex") rstyle(none)

end/*%>*/

capture program drop graph_p9099_wealth/*%<*/
program define graph_p9099_wealth

    /***************************************************************************
        Wealth concentration: 3 groups
    ***************************************************************************/

    /***********************************************************************
        Load wealth shares under capitalized specification and keep 
            groups we want; calculate wealth shares of each group; then save 
            as tempfiles.
    ***********************************************************************/

    foreach capitalized in preferred equreturns {

        if "`capitalized'" == "preferred" {
            load_analysis_data szz
        }
        else {
            load_analysis_data sz_v3
            rename hweal_szv3 hweal_equreturns
        }

        keep if inlist(group, "P0-90", "P90-99", "P99-100")

        isid group year

        #delimit ;
        replace group = cond(group == "P0-90", "_bot90", 
                        cond(group == "P90-99", "_p90_99", 
                        cond(group == "P99-100", "_top1", "alpaca")));
        #delimit cr
        assert group != "alpaca"

        sort year
        by year: egen tthweal_`capitalized' = total(hweal_`capitalized')

        gen wlthshare_`capitalized' = (hweal_`capitalized' / tthweal_`capitalized') * 100

        keep year group wlthshare_`capitalized'

        reshape wide wlthshare_`capitalized', i(year) j(group, string)

        tempfile `capitalized'
        save ``capitalized''
    }

    /***********************************************************************
        Load DFA data and calculate wealth shares
    ***********************************************************************/
            
    load_analysis_data dfaraw

    #delimit ;
    gen group = cond(inlist(category, "Bottom50", "Next40"), "_bot90", 
                cond(category == "Next9", "_p90_99", 
                cond(category == "Top1", "_top1", "alpaca")));
    #delimit cr
    assert group != "alpaca"    

    * Sum bottom two groups
    collapse (sum) networth, by(year date group)

    * Get annual averages from quarterly data
    collapse (mean) networth, by(year group)

    sort year
    by year: egen ttnetworth = total(networth)

    gen wlthshare_dfa = (networth / ttnetworth) * 100

    keep year group wlthshare_dfa

    reshape wide wlthshare_dfa, i(year) j(group, string)

    tempfile dfa
    save `dfa'

    /***********************************************************************
        Merge capitalized data together with DFA data
    ***********************************************************************/
        
    merge 1:1 year using `preferred'
    assert year > 2016 if _merge == 1
    assert year < 1989 if _merge == 2
    assert _merge == 3 if inrange(year, 1989, 2016)
    drop _merge 

    merge 1:1 year using `equreturns', assert(1 3)
    assert year > 2016 if _merge == 1
    assert _merge == 3 if inrange(year, 1966, 2016)
    drop _merge

    /***********************************************************************
        Plot series
    ***********************************************************************/

    colorpalette "$u3", intensity(0.1(.05)1)
    local u3l "`r(p6)'"
    colorpalette "$u1", intensity(0.1(.05)1)
    local u1l "`r(p6)'"
    colorpalette "$p2", intensity(0.1(.05)1)
    local p2l "`r(p6)'"

    keep if inrange(year, 1966, 2016)

    sort year
    #delimit ;
    graph twoway 
        (connected wlthshare_equreturns_bot90 year, msymbol(t) mc("`p2l'") lc("`p2l'"))
        (connected wlthshare_equreturns_p90_99 year, msymbol(dh) mc("`u3l'") lc("`u3l'"))
        (connected wlthshare_equreturns_top1 year, msymbol(o) mc("`u1l'") lc("`u1l'"))
        (connected wlthshare_preferred_bot90 year, msymbol(t) mc("$p2") lc("$p2"))
        (connected wlthshare_preferred_p90_99 year, msymbol(dh) mc("$u3") lc("$u3"))
        (connected wlthshare_preferred_top1 year, msymbol(o) mc("$u1") lc("$u1"))
        ,
        $gpr
        ytitle("Wealth Share (%)") xtitle("")
        xlab(1965(5)2015) xscale(range(1964 2016))
        yscale(range(17.5 47.5))
        legend(order(
            4 "P0-90 Baseline"      5 "P90-99 Baseline"  
            6 "P99-100 Baseline"    1 "P0-90 Equal Returns"  
            2 "P90-99 Equal Returns"      3 "P99-100 Equal Returns") 
            region(lcolor(white) margin(tiny)) rows(2) symxsize(*.5) size(medsmall))
        xsize(7.5);
    #delimit cr
    graph_export "wealth_comparison_9099"

end/*%>*/

capture program drop graph_topwealth_cross /*%<*/
program define graph_topwealth_cross

    /***************************************************************************
        Prepare SCF + Forbes 400 data 
    ***************************************************************************/

    /***********************************************************************
        Load SCF bulletin file for 2016 keep risk-based portfolio
            aggregates, and get totals by ES-adjusted within-decile groups
    ***********************************************************************/
    load_analysis_data scfcross_v3
    drop if grp == 0
    tempfile scf
    save `scf'

    /***********************************************************************
        Prepare Forbes wealth using Laurence classification of 
            public and private equity wealth
    ***********************************************************************/ 
    load_analysis_data scfcross_v3_w400
    drop if grp == 0
    tempfile scf_f400
    save `scf_f400'

    /***************************************************************************
        Prepare tax data shares of aggregates
    ***************************************************************************/

    /***********************************************************************
        Use program to load aggregate wealth under preferred and 
            equal-returns specification
    ***********************************************************************/
    load_analysis_data szz

    keep if year == 2016 & group == "All"
    keep year hweal$sz_wlth_defn hweal$sz_wlth_defn_late_v3 hweal$preferred_defn_late
    assert _N == 1

    global hweal_pref = hweal$preferred_defn_late[1]
    global hweal_equrtrns = hweal$sz_wlth_defn[1]
    global hweal_equ = hweal$sz_wlth_defn_late_v3[1]
    global hweal_base = hweal$preferred_defn_late[1]
    
    /***********************************************************************
        Load wealth by portfolio category 
    ***********************************************************************/
    load_analysis_data telescope_parts_v3

    *foreach portfoliocat in hwbus hwequ taxbond pen_hou_oth {
    foreach portfoliocat in hwfix hwbus pthru hwequ ccorw hwpen hwhou hwoth pen_hou_oth {
        gen `portfoliocat'_prefAsh = (`portfoliocat'_pref / $hweal_pref) * 100
        gen `portfoliocat'_equrtrnsAsh = (`portfoliocat'_equrtrns / $hweal_equrtrns) * 100
        gen `portfoliocat'_baseAsh = (`portfoliocat'_base / $hweal_base) * 100
        gen `portfoliocat'_equAsh = (`portfoliocat'_equ / $hweal_equ) * 100
    }

    keep grp *Ash

    tempfile inside_collapse
    save `inside_collapse'

    /***************************************************************************
        Get DFA data
    ***************************************************************************/
    load_analysis_data dfacross
    tempfile dfa
    save `dfa'

    /***************************************************************************
        Merge SCF together with tax data and DFA; make graphs
    ***************************************************************************/
    use `scf_f400', clear
    merge 1:1 grp using `inside_collapse', keep(3) nogen
    merge 1:m grp using `dfa', keep(3) nogen

    sort grp

    export delimited using wealth_composition_telescope.csv, replace

    foreach portfoliocat in hwfix hwbus pthru hwequ ccorw hwpen hwhou hwoth pen_hou_oth {

        #delimit ;
        twoway (connect `portfoliocat'_equAsh grp, ms(s) lc("$u3") mc("$u3"))
            (connect `portfoliocat'_baseAsh grp, ms(O) lc("$u1") mc("$u1") lwidth(thick))
            (connect `portfoliocat'_scfAsh grp, ms(Th) lc("$p2") mc("$p2") lw(thin))
            (connect `portfoliocat'_scf_f400Ash grp, 
                ms(T) lc("$p2") mc("$p2") lw(medthin))
            (connect `portfoliocat'_dfaAsh grp if category == "Next9", 
                ms(none) lc(black) lpattern(-))
            (connect `portfoliocat'_dfaAsh grp if category == "Top1", 
                ms(none) lc(black) lpattern(_))
            ,
            xlab(1 "P90-91" 2 "P91-92" 3 "P92-93" 4 "P93-94" 5 "P94-95" 
                 6 "P95-96" 7 "P96-97" 8 "P97-98" 9 "P98-99" 10 "P99-99.9" 
                 /* 11 "P99.9-100" */ 
                 11 "P99.9-99.99" 12 "P99.99-100", angle(30) labsize(small))
            /* yscale(range(0 11)) ylab(0(2)10) */
            yscale(range(0 8)) ylab(0(2)8) 
            xtitle("")
            ytitle("Share of Total Household Wealth (%)") $gpr
            legend(order(
                         4 "Harmonized SCF w/Forbes"
                         1 "Equal Returns"
                         5 "DFA P90-99 (Even)"
                         3 "Harmonized SCF"
                         2 "Baseline"
                         6 "DFA Top 1% (Even)") size(small) symxsize(*.5)
                    col(3) region(lc(white) margin(tiny)))
            xsize(6);
        #delimit cr
        graph_export `portfoliocat'_share_agg_wealth
    }

    #delimit ;
    twoway (connect ccorw_equAsh grp, ms(s) color("$u3"))
        (connect ccorw_baseAsh grp, ms(O) color("$u1") lwidth(thick))
        (connect ccorw_scfAsh grp, ms(Th) color("$p2") lwidth(thin))
        (connect privccorw_scfAsh grp, ms(Th) color("$f1") lwidth(thin) lp("-"))
        (connect ccorw_scf_f400Ash grp, 
            ms(T) color("$p2") lwidth(medthin) lp(_))
        (connect ccorw_dfaAsh grp if category == "Next9", 
            ms(none) lc(black) lpattern(-))
        (connect ccorw_dfaAsh grp if category == "Top1", 
            ms(none) lc(black) lpattern(_))
        ,
        xlab(1 "P90-91" 2 "P91-92" 3 "P92-93" 4 "P93-94" 5 "P94-95" 
             6 "P95-96" 7 "P96-97" 8 "P97-98" 9 "P98-99" 10 "P99-99.9" 
             11 "P99.9-P99.99" 12 "P99.99-100", angle(30) labsize(small))
        yscale(range(0 11)) ylab(0(2)10)
        xtitle("")
        ytitle("Share of Total Household Wealth (%)") $gpr
        legend(order(5 "Harmonized SCF w/Forbes"
                     1 "Equal Returns"
                     6 "DFA P90-99 (Even)"
                     3 "Harmonized SCF"
                     2 "Baseline"
                     7 "DFA Top 1% (Even)"
                     4 "Harmonized SCF: Private C-corps only") 
                     size(small) symxsize(*.5)
                col(3) region(lc(white) margin(tiny)))
        xsize(6);
    #delimit cr
    graph_export ccorw_priv_breakout_share_agg_wealth

end/*%>*/

capture program drop graph_components/*%<*/
program define graph_components

    /***************************************************************************
        Aggregate Top 0.1% Wealth Composition
    ***************************************************************************/

    /***********************************************************************
        Get shares from family offices (source = UBS Global Family 
            Office Report for 2016, saved in ~/Dropbox/wealth/data/
            PIMCO/UBS family office reports) [this is in progress for now]
    ***********************************************************************/

    load_analysis_data family
    keep hw* ccorw_sh pthru_sh pen_oth_sh

    gen dataspec = "ubsfamilyoffice"

    tempfile ubsfamilyoffice
    save `ubsfamilyoffice'

    /***********************************************************************
        Prepare estate tax data (these are already in shares)
    ***********************************************************************/

    load_analysis_data estatecross_mgs
    keep *_sh

    gen pen_oth_sh = hwpen_sh + hwoth_sh

    foreach sharevar of varlist _all {
        replace `sharevar' = 100 * `sharevar'
    }

    gen dataspec = "estate"

    tempfile estate_portfshares
    save `estate_portfshares'

    /***********************************************************************
        Retrieve preferred capitalized wealth of top 0.1% with 
            five portfolio categories
    ***********************************************************************/

    foreach wlthgrp in top10 top1 top01 top001 top0001 {

        load_analysis_data szz 
        keep if year == 2016

        #delimit ;
        local w_group = cond("`wlthgrp'" == "top10", 2, 
                        cond("`wlthgrp'" == "top1", 3, 
                        cond("`wlthgrp'" == "top01", 4, 
                        cond("`wlthgrp'" == "top001", 5, 6))));
        #delimit cr

        keep if w${preferred_defn_late}_group == `w_group'

        assert _N == 1

        gen pen_oth = hwpen_base + hwoth_base

        egen hweal_check1 = rowtotal(hwfix_base hwbus_base hwhou_base pen_oth hwequ_base)
        egen hweal_check2 = rowtotal(hwfix_base pthru_base hwhou_base pen_oth ccorw_base)
        
        if "`wlthgrp'" != "top0001" {
            assert inrange(hweal_check1 / hweal$preferred_defn_late, 0.99, 1.01) 
            assert inrange(hweal_check2 / hweal$preferred_defn_late, 0.99, 1.01) 
        }
        else { // These have to be censored for privacy reasons
            assert inrange(hweal_check1 / hweal$preferred_defn_late, 0.95, 1.05)  
            assert inrange(hweal_check2 / hweal$preferred_defn_late, 0.95, 1.05)  

            * Due to rounding, going to use the components for denominator
            replace hweal$preferred_defn_late = hweal_check2
        }
        keep hweal$preferred_defn_late hw*_base pthru_base ccorw_base pen_oth

        rename (hweal$preferred_defn_late *_base) (hweal *)

        gen dataspec = "preferred`wlthgrp'"

        tempfile preferred`wlthgrp'
        save `preferred`wlthgrp''
    }

    /***********************************************************************
        Retrieve equal returns capitalized wealth of top groups 
            with five portfolio categories
    ***********************************************************************/

    foreach wlthgrp in top10 top1 top01 top001 top0001 {

        load_analysis_data sz_v3
        keep if year == 2016

        #delimit ;
        local w_group = cond("`wlthgrp'" == "top10", 2, 
                        cond("`wlthgrp'" == "top1", 3, 
                        cond("`wlthgrp'" == "top01", 4, 
                        cond("`wlthgrp'" == "top001", 5, 6))));
        #delimit cr

        keep if w${sz_wlth_defn_late_v3}_group == `w_group'

        assert _N == 1

        gen pen_oth_equ = hwpen_equ + hwoth_equ
        
        qui egen hweal_chk1 = rowtotal(hwfix_equ hwbus_equ hwhou_equ pen_oth_equ hwequ_equ)
        qui egen hweal_chk2 = rowtotal(hwfix_equ pthru_equ hwhou_equ pen_oth_equ ccorw_equ)

        #delimit ;
        qui egen hweal_chk3 = rsum(ccorw scorw taxbond taxbond_mufmisc_sz muni currency 
                                   rentalhome ownerhome partw_sz20_scaled
                                   solepropw_sz20_scaled 
                                   hwpen_ini ownermort rentalmort nonmort_ini);
        #delimit cr

        if "`wlthgrp'" != "top0001" { 
            assert inrange(hweal_chk1 / hweal$sz_wlth_defn_late_v3, 0.9999, 1.0001)
            assert inrange(hweal_chk2 / hweal$sz_wlth_defn_late_v3, 0.9999, 1.0001)
            assert inrange(hweal_chk3 / hweal$sz_wlth_defn_late_v3, 0.9999, 1.0001)
        }
        else { // These have to be censored for privacy reasons
            assert inrange(hweal_chk1 / hweal$sz_wlth_defn_late_v3, 0.925, 1.075)
            assert inrange(hweal_chk2 / hweal$sz_wlth_defn_late_v3, 0.925, 1.075)
            assert inrange(hweal_chk3 / hweal$sz_wlth_defn_late_v3, 0.925, 1.075)

        }

        keep hweal$sz_wlth_defn_late_v3 hwfix_equ hwequ_equ hwbus_equ hwpen_equ ///
            hwhou_equ hwoth_equ pthru_equ ccorw_equ pen_oth_equ

        rename *_equ *

        rename hweal$sz_wlth_defn_late_v3 hweal

        gen dataspec = "equreturns`wlthgrp'"

        tempfile equreturns`wlthgrp'
        save `equreturns`wlthgrp''
    }

    /***********************************************************************
        Get top group totals from SCF
    ***********************************************************************/
    load_analysis_data scfrevision_v3

    gen privcorp = privccorw + scorw

    assert othdebt <= 0
    replace hwoth = hwoth + othdebt
    gen pen_oth = hwpen + hwoth

    ds hwequ hwfix hwbus hwpen hwhou hwoth pen_oth ccorw pthru privccorw privcorp
    local othersplitvars "`r(varlist)'"

    es_rank_scf, rankvar(networth_pref) othersplitvars("`othersplitvars'")

    tempfile orig 
    save `orig'

    foreach wlthgrp in top1 top01 top001 top0001 {

        use `orig', clear

        #delimit ;
        local threshold = cond("`wlthgrp'" == "top1", 0.99, 
                          cond("`wlthgrp'" == "top01", 0.999, 
                          cond("`wlthgrp'" == "top001", 0.9999, 0.99999)));
        #delimit cr

        keep if es_rank >= `threshold' & year == 2016

        egen networth_check1 = rowtotal(hwequ hwfix hwbus hwhou pen_oth)
        egen networth_check2 = rowtotal(ccorw hwfix pthru hwhou pen_oth)
        assert inrange(networth_pref / networth_check1, 0.999, 1.001) & ///
            inrange(networth_pref / networth_check2, 0.999, 1.001)

        collapse (sum) hweal = networth_pref `othersplitvars' [fw = wgt1B]

        assert _N == 1

        foreach billionX of varlist hweal `othersplitvars' {
            replace `billionX' = `billionX' / 1E9
        }

        if "`wlthgrp'" == "top001" {
            foreach sharevar of varlist hwbus hwequ hwfix hwhou pen_oth {
                local `sharevar'_sh_`wlthgrp' = (`sharevar' / hweal)
                assert inrange(``sharevar'_sh_`wlthgrp'', 0, 1)
            }

            local privccorw_sh_hwbus = privccorw / hwbus 
            local ccorw_sh_privcorp = privccorw / privcorp 

            assert inrange(`privccorw_sh_hwbus', 0, 1) & inrange(`ccorw_sh_privcorp', 0, 1)
        }
        drop privccorw

        gen dataspec = "scf`wlthgrp'"

        tempfile scf`wlthgrp'
        save `scf`wlthgrp''
    }

    /***************************************************************************
        Get top 1% totals from DFA
    ***************************************************************************/
    use $inputs/dfa_revision.dta, clear

    keep if year == 2016 & category == "Top1"
    assert _N == 1
    drop year category

    assert sign(othdebt) == -1
    replace hwoth_dfa = hwoth_dfa + othdebt_dfa
    gen pen_oth = hwpen_dfa + hwoth_dfa

    keep networth_pref hwfix_dfa hwbus_dfa pthru_dfa hwequ_dfa ccorw_dfa ///
        hwhou_dfa hwpen_dfa hwoth_dfa pen_oth
    rename (*_dfa networth_pref) (* hweal)

    foreach unscaled of varlist _all { // Scale out of millions into dollars
        replace `unscaled' = `unscaled' * 1E6
    }

    gen dataspec = "dfa"

    tempfile dfa
    save `dfa'


    /***********************************************************************
        Allocate Forbes total wealth according to top 0.01% in 
            SCF
    ***********************************************************************/

    import delim forbesportfolios.csv, clear

    gen pen_oth = hwpen_f400 + hwoth_f400
    drop pen_hou_oth_f400

    rename *_f400 *

    gen dataspec = "forbes"

    tempfile forbes
    save `forbes'

    /***********************************************************************
        Merge together capitalized and SCF; calculate portfolio 
            shares using these aggregates
    ***********************************************************************/

    use `preferredtop1', clear

    #delimit ;
    append using `preferredtop01' `preferredtop001' `preferredtop0001' 
        `equreturnstop1' `equreturnstop01' `equreturnstop001' `equreturnstop0001' 
        `scftop1' `scftop01' `scftop001' `scftop0001' `forbes' `dfa';
    #delimit cr

    keep hw* ccorw pthru pen_oth dataspec

    local lastrow = _N
    forv row = 1 / `lastrow' {
        local newname`row' = dataspec[`row']
    }

    xpose, clear varname
    qui drop if _varname == "dataspec"
    forv column = 1 / `lastrow' {
        rename v`column' `newname`column''
    }

    qui gen scftop1plusforbes = scftop1 + forbes
    qui gen scftop01plusforbes = scftop01 + forbes
    qui gen scftop001plusforbes = scftop001 + forbes
    qui gen scftop0001plusforbes = scftop0001 + forbes
    drop forbes

    local lastrow = _N
    forv row = 1 / `lastrow' {
        local v`row' = _varname[`row']
    }

    qui export delim figure13xpose.csv, replace // for portfolio totals appx fig

    drop _varname

    xpose, clear varname

    forv column = 1 / `lastrow' {
        local newname = "`v`column''"
        rename v`column' `newname'
    }
    rename _varname dataspec
    drop hwpen hwoth


    foreach sharevar of varlist hwequ ccorw hwfix hwbus pthru hwhou pen_oth {
        qui gen `sharevar'_sh = (`sharevar' / hweal) * 100
    }

        /***********************************************************************
            Append shares from other sources, where portfolio shares 
                are pre-calculated
        ***********************************************************************/

    append using `ubsfamilyoffice' `estate_portfshares'

        /***********************************************************************
            Make graphs
        ***********************************************************************/

    foreach grp in top1 top01 top001 top0001 {

        #delimit ;
        local lastsource = cond("`grp'" == "top1", "dfa",
                           cond("`grp'" == "top01", "estate", "ubsfamilyoffice"));

        local lastsrclab = cond("`grp'" == "top1", "DFA",
                           cond("`grp'" == "top01", "Estate Tax", "UBS Fmly Offc"));

        qui gen `grp'_barorder = cond(dataspec == "equreturns`grp'", 1, 
                                 cond(dataspec == "preferred`grp'", 2, 
                                 cond(dataspec == "scf`grp'plusforbes", 3, 
                                 cond(dataspec == "`lastsource'", 4, .))));
        #delimit cr

        #delimit ;
        graph hbar hwfix_sh hwequ_sh hwbus_sh hwhou_sh pen_oth_sh
            if !missing(`grp'_barorder)
            ,
            $gpr stack 
            over(`grp'_barorder,
                relabel(1 "Equal Returns" 2 "Baseline" 3 "SCF + Forbes" 4 "`lastsrclab'") 
                lab(labsize(small)))
            legend(order(1 "Fixed Income" 2 "Public Equity" 3 "Private Business"
                         4 "Housing" 5 "Pensions and Other")
                    size(small) region(lcolor(white) margin(tiny)) row(2)
                    symxsize(*0.4))
            blabel(bar, position(center) color(white) format(%9.1fc) size(small))
            bar(1, fcolor("$u8") lcolor(black)) 
            bar(2, color("$u7")) 
            bar(3, color("$p21"))
            bar(4, color("$u9"))
            bar(5, color("$u10")) 
            xsize(6);
        graph_export `grp'portfolioshares;

        graph hbar hwfix_sh ccorw_sh pthru_sh hwhou_sh pen_oth_sh
            if !missing(`grp'_barorder)
            ,
            $gpr stack 
            over(`grp'_barorder,
                relabel(1 "Equal Returns" 2 "Baseline" 3 "SCF + Forbes" 4 "`lastsrclab'")
                lab(labsize(small)))
            legend(order(1 "Fixed Income" 2 "C-corp" 3 "Pass-through"
                         4 "Housing" 5 "Pensions, Other")
                    size(small) region(lcolor(white) margin(tiny)) row(1)
                    symxsize(*0.4))
            blabel(bar, position(center) color(white) format(%9.1fc) size(small))
            bar(1, fcolor("$u8") lcolor(black)) 
            bar(2, color("$u7")) 
            bar(3, color("$p21"))
            bar(4, color("$u9"))
            bar(5, color("$u10"))
            xsize(6);
        #delimit cr
        graph_export `grp'portfolioshares_ccorw_pthru
    }

end/*%>*/

capture program drop graph_components_overtime/*%<*/
program define graph_components_overtime
        
    load_analysis_data cpi
    tempfile cpi
    save `cpi'
    /***************************************************************************
        Retrieve top 1% and top 0.1% preferred capitalized portfolio 
            shares
    ***************************************************************************/
 
    load_analysis_data szz

    keep if inlist(w_group, 1, 7, 8, 3, 4, 5)

    assert inlist(group, "All", "P0-90", "P90-99", "P99-100", "P99.9-100", "P99.99-100")

    assert sign(nonmort) == -1
    gen pen_hou_oth_pref = hwpen_pref + hwhou_pref + nonmort + miscw
    gen pen_hou_oth_base = hwpen_base + hwhou_base + nonmort_ini + miscw

    gen hweal_prefH = hweal_preferred * (group == "All")
    bys year: egen total_hweal_preferred = max(hweal_prefH)
    drop hweal_prefH
    gen miss_prefH = missing_scorp * (group == "All")
    bys year: egen total_missing_scorp = max(miss_prefH)
    drop miss_prefH
    gen miss_prefH = missing_pship * (group == "All")
    bys year: egen total_missing_pship = max(miss_prefH)
    drop miss_prefH
    drop if group == "All"

    foreach sharevar of varlist *_pref *_base {
        gen `sharevar'_sh = (`sharevar' / hweal_preferred) * 100
        gen `sharevar'_tsh = (`sharevar' / total_hweal_preferred) * 100
    }

    keep year group *_sh *_tsh *_pref *_base

    rename (*_base_sh *_base_tsh *_base *_pref_sh *_pref_tsh *_pref) ///
         (*_sh_preferred *_tsh_preferred *_preferred *_sh_supple *_tsh_supple *_supple)

    replace group = cond(group == "P99-100", "_top1", ///
                    cond(group == "P99.9-100", "_top01", ///
                    cond(group == "P90-99", "_p9099", ///
                    cond(group == "P0-90", "_bot90", "_top001"))))

    reshape wide *_preferred *_supple, i(year) j(group, string)

    merge 1:1 year using `cpi', keep(3) keepusing(adjfactor2016) nogen
    foreach infladjvar of varlist hwfix_preferred* hwequ_preferred* hwbus_preferred* hwhou_preferred* pthru_preferred* ccorw_preferred* hwpen_preferred* pen_hou_oth_preferred* *_supple* {
        replace `infladjvar' = 1e-12 * `infladjvar' * adjfactor2016
    }

    tempfile preferred
    save `preferred'

    /***************************************************************************
        Retrieve top 1% and top 0.1% preferred capitalized portfolio 
            shares for equal returns with equal returns ranks
    ***************************************************************************/

    load_analysis_data sz_v3

    gen w_group = cond(year < 1980, w${sz_wlth_defn_early_v3}_group, ///
                        cond(year >= 1980 & year < 2001, w${sz_wlth_defn_mid_v3}_group, ///
                            w${sz_wlth_defn_late_v3}_group))

    keep if inlist(w_group, 1, 7, 8, 3, 4, 5)

    assert inlist(group, "All", "P0-90", "P90-99", "P99-100", "P99.9-100", "P99.99-100")

    assert sign(nonmort_ini) == -1
    gen pen_hou_oth_equ = hwpen_equ + hwhou_equ + nonmort_ini + miscw

    gen hweal_equH = hweal_szv3 * (group == "All")
    bys year: egen total_hweal_szv3 = max(hweal_equH)
    drop hweal_equH
    drop if group == "All"

    foreach sharevar of varlist *_equ {
        gen `sharevar'_sh = (`sharevar' / hweal_szv3) * 100
        gen `sharevar'_tsh = (`sharevar' / total_hweal_szv3) * 100
    }

    keep year group *_sh *_tsh *_equ 

    rename (*_equ_sh *_equ_tsh *_equ) (*_sh_equ *_tsh_equ *_equ)

    replace group = cond(group == "P99-100", "_top1", ///
                    cond(group == "P99.9-100", "_top01", ///
                    cond(group == "P90-99", "_p9099", ///
                    cond(group == "P0-90", "_bot90", "_top001"))))

    reshape wide *_equ, i(year) j(group, string)

    merge 1:1 year using `cpi', keep(3) keepusing(adjfactor2016) nogen
    foreach infladjvar of varlist hwfix_equ* hwhou_equ* pthru_equ* ccorw_equ* hwpen_equ* {
        replace `infladjvar' = 1e-12 * `infladjvar' * adjfactor2016
    }

    tempfile szv3
    save `szv3'

    /***************************************************************************
        Retrieve top 1% and top 0.1% preferred equal returns portfolio 
            shares
    ***************************************************************************/

    load_analysis_data sz

    keep if inlist(w${sz_wlth_defn}_group, 1, 7, 8, 3, 4, 5)
    assert inlist(group, "All", "P0-90", "P90-99", "P99-100", "P99.9-100", "P99.99-100")
    
    rename nonmort_ini hwoth
    gen pen_oth = hwpen_ini + hwoth

    gen hwequ = 0.8 * ccorw_ini
    gen hwbus = scorw_ini + partw_ini + solepropw_ini + 0.2 * ccorw_ini
    gen hwfix = taxbond_ini + muni_ini + currency_ini
    gen pthru = scorw_ini + partw_ini + solepropw_ini
    drop ccorw
    gen ccorw = ccorw_ini
    
    drop hwhou 
    gen hwhou = rentalhome_ini + ownerhome_ini + ownermort_ini + rentalmort_ini

    gen pen_hou_oth = hwpen_ini + hwhou + hwoth

    drop hwpen
    gen hwpen = hwpen_ini

    egen hweal_check = rowtotal(hwfix hwbus hwequ pen_hou_oth)
    assert inrange(hweal_check / hweal$sz_wlth_defn, 0.999, 1.001)

    gen hweal${sz_wlth_defn}H = hweal${sz_wlth_defn} * (group == "All")
    bys year: egen total_hweal${sz_wlth_defn} = max(hweal${sz_wlth_defn}H)
    drop hweal${sz_wlth_defn}H
    drop if group == "All"

    foreach sharevar in hwfix hwbus hwequ hwpen hwhou hwoth pthru ccorw pen_hou_oth {
        gen `sharevar'_sh = (`sharevar' / hweal$sz_wlth_defn) * 100
        gen `sharevar'_tsh = (`sharevar' / total_hweal$sz_wlth_defn) * 100
    }

    keep year group *_sh *_tsh hwfix hwbus hwequ hwpen hwhou hwoth pthru ccorw pen_hou_oth

    rename (*_sh *_tsh hwfix hwbus hwequ hwpen hwhou pthru ccorw hwoth pen_hou_oth) ///
        (*_sh_equal *_tsh_equal hwfix_equal hwbus_equal hwequ_equal hwpen_equal hwhou_equal pthru_equal ccorw_equal hwoth_equal pen_hou_oth_equal)

    replace group = cond(group == "P99-100", "_top1", ///
                    cond(group == "P99.9-100", "_top01", ///
                    cond(group == "P90-99", "_p9099", ///
                    cond(group == "P0-90", "_bot90", "_top001"))))

    reshape wide *_equal, i(year) j(group, string)

    merge 1:1 year using `cpi', keep(3) keepusing(adjfactor2016) nogen

    foreach infladjvar of varlist hwfix_equal* hwequ_equal* hwbus_equal* hwhou_equal* hwpen_equal* pthru_equal* ccorw_equal* hwoth_equal* pen_hou_oth_equal* {
        replace `infladjvar' = 1e-12 * `infladjvar' * adjfactor2016
    }

    tempfile equreturns
    save `equreturns'

    /***************************************************************************
        Retrieve SCF portfolio shares
    ***************************************************************************/

        /***********************************************************************
            Get 2016 Forbes public-private, C corp-pass through 
                allocations; these are based on individual Forbes hunting done
                by Laurence O'Brien, so we think they're better than just
                allocating according to tippy-top SCF work
        ***********************************************************************/

    import delimited forbesportfolios.csv, clear    

    gen f400biz = ccorw_f400 + pthru_f400
    gen f400biz_check = hwequ_f400 + hwbus_f400
    assert inrange(f400biz / f400biz_check, 0.99999, 1.00001)

    local ccorw_sh_f400biz = ccorw_f400 / f400biz
    local pthru_sh_f400biz = 1 - `ccorw_sh_f400biz'

    local hwequ_sh_f400biz = hwequ_f400 / f400biz
    local hwbus_sh_f400biz = 1 - `hwequ_sh_f400biz'

        /***********************************************************************
            Get aggregate wealth by wealth group, adjusting ranks to 
                be compatible with equal-split counts
        ***********************************************************************/

    load_analysis_data scfrevision_v3

    es_rank_scf, rankvar(networth_pref) ///
        othersplitvars("hwequ ccorw hwfix hwbus pthru hwpen hwhou hwoth othdebt")

    /* Wealth groups: 6 = below P90, 5 = P90-99; 4 = P99-99.9; 3 = P99.9-99.99, 
        2 = P99.99-99.999, 1 = top 0.001% (P99.999-100) */
    #delimit ;
    gen wlthgrp = cond(es_rank < 0.90, 6, 
                  cond(es_rank < 0.99, 5, 
                  cond(es_rank < 0.999, 4, 
                  cond(es_rank < 0.9999, 3, 
                  cond(es_rank < 0.99999, 2, 1)))));
    #delimit cr

    assert sign(othdebt) == -1 | othdebt == 0
    replace hwoth = hwoth + othdebt

    #delimit ;
    collapse (sum) hweal_scf = networth_pref hwequ_scf = hwequ ccorw_scf = ccorw 
        hwfix_scf = hwfix hwbus_scf = hwbus pthru_scf = pthru hwpen_scf = hwpen 
        hwhou_scf = hwhou hwoth_scf = hwoth [fw = wgt1B], by(year wlthgrp);
    #delimit cr

    foreach billionX of varlist *_scf {
        replace `billionX' = `billionX' / 1E9 // Unscale from 1B weights
    }

    gen pen_hou_oth_scf = hwpen_scf + hwhou_scf + hwoth_scf

    gen hweal_check1 = hwequ_scf + hwfix_scf + hwbus_scf + pen_hou_oth_scf
    gen hweal_check2 = ccorw_scf + hwfix_scf + pthru_scf + pen_hou_oth_scf
    assert inrange(hweal_check1 / hweal_scf, 0.999, 1.001) & ///
           inrange(hweal_check2 / hweal_scf, 0.999, 1.001)
    drop hweal_check?

    reshape wide *_scf, i(year) j(wlthgrp)

    foreach wlthcmpt in hweal_scf hwequ_scf ccorw_scf hwfix_scf hwbus_scf ///
        pthru_scf hwpen_scf hwhou_scf hwoth_scf pen_hou_oth_scf {

        egen `wlthcmpt'_total = rowtotal(`wlthcmpt'?)

        gen `wlthcmpt'_top001 = `wlthcmpt'1 + `wlthcmpt'2
        gen `wlthcmpt'_top01 = `wlthcmpt'_top001 + `wlthcmpt'3
        gen `wlthcmpt'_top1 = `wlthcmpt'_top01 + `wlthcmpt'4
        
        drop `wlthcmpt'2 `wlthcmpt'3 `wlthcmpt'4
        
        rename (`wlthcmpt'6 `wlthcmpt'5 `wlthcmpt'1) ///
            (`wlthcmpt'_bot90 `wlthcmpt'_p9099 `wlthcmpt'_t0001)
    }

        /***********************************************************************
            Calculate top 0.01% portfolio shares for non-business 
                variables; will use these to allocate Forbes wealth across 
                categories
        ***********************************************************************/

    foreach nonbizcat in hwfix_scf hwpen_scf hwhou_scf hwoth_scf pen_hou_oth_scf {
        gen `nonbizcat'_top001share = `nonbizcat'_top001 / hweal_scf_top001
        
        if "`nonbizcat'" != "hwoth_scf" { // mainly non-mortgage debt, so negative
            assert inrange(`nonbizcat'_top001share, 0, 1)
        }
        else {
            assert inrange(abs(`nonbizcat'_top001share), 0, 1)
        }
    }

    gen biz_scf_top001share = (hwbus_scf_top001 + hwequ_scf_top001) / hweal_scf_top001
    gen biz_scf_top001share_check = (pthru_scf_top001 + ccorw_scf_top001) / hweal_scf_top001

    assert abs(biz_scf_top001share - biz_scf_top001share_check) < 1E-7
    drop biz_scf_top001share_check

        /***********************************************************************
            Merge on Forbes wealth and allocate across portfolio 
                categories
        ***********************************************************************/

    merge 1:1 year using $dtadir/forbeswlth.dta, assert(3) nogen

    rename forbeswlth hweal_f400

    foreach wlthcmpt in biz hwfix pen_hou_oth hwpen hwhou hwoth {
        gen `wlthcmpt'_f400 = `wlthcmpt'_scf_top001share * hweal_f400
        
        assert !missing(`wlthcmpt'_f400)
        drop `wlthcmpt'_scf_top001share
    }

    foreach bizcat in hwbus hwequ ccorw pthru {
        gen `bizcat'_f400 = ``bizcat'_sh_f400biz' * biz_f400
        assert !missing(`bizcat'_f400)
    }

    drop biz_f400

        /***********************************************************************
            Add Forbes wealth to appropriate SCF aggregates
        ***********************************************************************/

    foreach addf400 of varlist *_total *_top1 *_top01 *_top001 *_t0001 {
        local wlthcmpt = substr("`addf400'", 1, strpos("`addf400'", "scf") - 2)

        qui replace `addf400' = `addf400' + `wlthcmpt'_f400
    }

        /***********************************************************************
            Calculate wealth shares, factoring in Forbes
        ***********************************************************************/

    foreach wlthcmpt in hwequ ccorw hwfix hwbus pthru pen_hou_oth hwpen hwhou hwoth {
        foreach grp in bot90 p9099 top1 top01 top001 t0001 {
            gen `wlthcmpt'_sh_scf_`grp' = (`wlthcmpt'_scf_`grp' / hweal_scf_`grp') * 100
            gen `wlthcmpt'_tsh_scf_`grp' = (`wlthcmpt'_scf_`grp' / hweal_scf_total) * 100
            
            replace `wlthcmpt'_scf_`grp' = `wlthcmpt'_scf_`grp' / 1E12
        }    
    }

    tempfile scf
    save `scf'

    /***************************************************************************
        Retrieve DFA portfolio shares (top 1% only)
    ***************************************************************************/

    load_analysis_data dfarevision

    assert sign(othdebt_dfa) == -1
    qui replace hwoth_dfa = hwoth_dfa + othdebt_dfa

    qui gen pen_hou_oth_dfa = hwpen_dfa + hwhou_dfa + hwoth_dfa

    isid year category
    levelsof category, local(categories) clean
    assert "`categories'" == "Bottom50 Next40 Next9 Top1"

    gen group = cond(inlist(category, "Bottom50", "Next40"), "_bot90", ///
                cond(category == "Next9", "_p9099", "_top1"))

    collapse (sum) networth_pref hwequ_dfa ccorw_dfa hwfix_dfa hwbus_dfa pthru_dfa ///
        hwhou_dfa hwpen_dfa hwoth_dfa pen_hou_oth_dfa, by(year group)

    sort year group
    by year: egen total_wealth = total(networth_pref)

    qui gen hweal_check1 = hwequ_dfa + hwfix_dfa + hwbus_dfa + pen_hou_oth_dfa 
    qui gen hweal_check2 = ccorw_dfa + hwfix_dfa + pthru_dfa + pen_hou_oth_dfa 
    assert inrange(hweal_check1 / networth_pref, 0.999, 1.001) & ///
           inrange(hweal_check2 / networth_pref, 0.999, 1.001)

    foreach wlthcmpt in hwequ ccorw hwfix hwbus pthru hwhou hwpen hwoth pen_hou_oth {
        qui gen `wlthcmpt'_sh_dfa = (`wlthcmpt'_dfa / networth_pref) * 100
        qui gen `wlthcmpt'_tsh_dfa = (`wlthcmpt'_dfa / total_wealth) * 100

        replace `wlthcmpt'_dfa = `wlthcmpt'_dfa / 1E6 // Scale out of millions into trillions
    }

    drop total_wealth networth_pref hweal_check?

    reshape wide *_dfa, i(year) j(group, string)

    tempfile dfa
    save `dfa'

    /***************************************************************************
        Merge together data from various sources
    ***************************************************************************/
    use `preferred', clear
    merge 1:1 year using `equreturns', assert(3) nogen
    merge 1:1 year using `szv3', assert(3) nogen
    
    merge 1:1 year using `scf'
    assert _merge == 3 if inrange(year, 1989, 2016) & mod(year - 1989, 3) == 0
    assert year == 2019 if _merge == 2
    assert !inrange(year, 1989, 2019) | mod(year - 1989, 3) > 0 if _merge == 1
    drop _merge

    merge 1:1 year using `dfa'
    assert _merge == 3 if inrange(year, 1989, 2016) | year == 2019
    assert inlist(year, 2017, 2018) if _merge == 2
    assert year < 1989 if _merge == 1
    drop _merge 

    outsheet using wealthcompovertime.csv, replace

    /***************************************************************************
        Make plots
    ***************************************************************************/
	
    sort year
    foreach sharevar in hwfix ccorw pthru hwhou hwpen {/*%<*/
    
        #delimit ;
        local cmptname = cond("`sharevar'" == "hwfix", "Fixed Income", 
                         cond("`sharevar'" == "hwequ", "Public Equity", 
                         cond("`sharevar'" == "hwbus", "Private Business", 
                         cond("`sharevar'" == "pthru", "Pass-Through Business", 
                         cond("`sharevar'" == "ccorw", "C-corporation Equity", 
                         cond("`sharevar'" == "hwpen", "Pensions", 
                         cond("`sharevar'" == "hwhou", "Housing", 
                                "Pensions, Housing, Other")))))));

        twoway (connect `sharevar'_sh_preferred_top1 year, ms(oh) lc("$u1") mc("$u1") lwidth(medthick))
            (connect `sharevar'_sh_equ_top1 year, ms(s) lc("$u3") mc("$u3")) 
            (connect `sharevar'_sh_scf_top1 year, ms(t) lc("$p2") mc("$p2") lwidth(thin))
            (connect `sharevar'_sh_dfa_top1 year, ms(d) lc("$u4") mc("$u4") lpattern("-."))
                if inrange(year, 1960, 2019)
            ,
            $gpr ytitle("`cmptname' Sh Top 1% Wealth (%)") xtitle("")
            legend(order(2 "Equal" 1 "Baseline" 3 "SCF" 4 "DFA") 
                    region(lcolor(white) margin(tiny)) row(1))
            ysca(range(0 70)) ylab(0(10)70)
            xsize(5.5);
        graph_export `sharevar'_shares_top1;

        twoway (connect `sharevar'_sh_preferred_top01 year, ms(oh) lc("$u1") mc("$u1") lwidth(medthick))
            (connect `sharevar'_sh_equ_top01 year,  ms(s) lc("$u3") mc("$u3")) 
            (connect `sharevar'_sh_scf_top01 year,  ms(t) lc("$p2") mc("$p2") lwidth(thin))
                if inrange(year, 1960, 2019)
            ,
            $gpr ytitle("`cmptname' Sh Top 0.1% Wealth (%)") xtitle("")
            legend(order(2 "Equal" 1 "Baseline" 3 "SCF") 
                    region(lcolor(white) margin(tiny)) row(1))
            ysca(range(0 70)) ylab(0(10)70)
            xsize(5.5);
        graph_export `sharevar'_shares_top01;

        twoway (connect `sharevar'_sh_preferred_top001 year, ms(oh) lc("$u1") mc("$u1") lwidth(medthick))
            (connect `sharevar'_sh_equ_top001 year,  ms(s) lc("$u3") mc("$u3")) 
            (connect `sharevar'_sh_scf_top001 year, ms(t) lc("$p2") mc("$p2") lwidth(thin))
                if inrange(year, 1960, 2019)
            ,
            $gpr ytitle("`cmptname' Sh Top 0.01% Wealth (%)") xtitle("")
            legend(order(2 "Equal" 1 "Baseline" 3 "SCF") 
                    region(lcolor(white) margin(tiny)) row(1))
            ysca(range(0 70)) ylab(0(10)70)
            xsize(5.5);
        graph_export `sharevar'_shares_top001;

        twoway (connect `sharevar'_sh_preferred_p9099 year, ms(oh) lc("$u1") mc("$u1") lwidth(medthick))
            (connect `sharevar'_sh_equ_p9099 year,  ms(s) lc("$u3") mc("$u3"))
            (connect `sharevar'_sh_scf_p9099 year, ms(t) lc("$p2") mc("$p2") lwidth(thin))
            (connect `sharevar'_sh_dfa_p9099 year, ms(d) lc("$u4") mc("$u4") lpattern("-."))
                if inrange(year, 1960, 2019)
            ,
            $gpr ytitle("`cmptname' Sh P90-99 Wealth (%)") xtitle("")
            legend(order(2 "Equal" 1 "Baseline" 3 "SCF" 4 "DFA") 
                    region(lcolor(white) margin(tiny)) row(1))
            ysca(range(0 70)) ylab(0(10)70)
            xsize(5.5);
        graph_export `sharevar'_shares_p9099;

        twoway (connect `sharevar'_sh_preferred_bot90 year, ms(oh) lc("$u1") mc("$u1") lwidth(medthick))
            (connect `sharevar'_sh_equ_bot90 year,  ms(s) lc("$u3") mc("$u3"))
            (connect `sharevar'_sh_scf_bot90 year, ms(t) lc("$p2") mc("$p2") lwidth(thin))
            (connect `sharevar'_sh_dfa_bot90 year, ms(d) lc("$u4") mc("$u4") lpattern("-."))
                if inrange(year, 1960, 2019)
            ,
            $gpr ytitle("`cmptname' Sh P0-90 Wealth (%)") xtitle("")
            legend(order(2 "Equal" 1 "Baseline" 3 "SCF" 4 "DFA") 
                    region(lcolor(white) margin(tiny)) row(1))
            ysca(range(0 70)) ylab(0(10)70)
            xsize(5.5);
        graph_export `sharevar'_shares_bot90;
        #delimit cr
    }

    foreach sharevar in hwfix ccorw pthru hwpen hwhou {
    
        #delimit ;
        local cmptname = cond("`sharevar'" == "hwfix", "Fixed Income", 
                         cond("`sharevar'" == "hwequ", "Public Equity", 
                         cond("`sharevar'" == "hwbus", "Private Business", 
                         cond("`sharevar'" == "pthru", "Pass-Through Business", 
                         cond("`sharevar'" == "ccorw", "C-corporation Equity", 
                         cond("`sharevar'" == "hwpen", "Pensions", 
                         cond("`sharevar'" == "hwhou", "Housing", 
                                "Pensions, Housing, Other")))))));

        twoway (connect `sharevar'_tsh_preferred_top1 year, ms(oh) lc("$u1") mc("$u1") lwidth(medthick))
            (connect `sharevar'_tsh_equ_top1 year, ms(s) lc("$u3") mc("$u3")) 
            (connect `sharevar'_tsh_scf_top1 year, ms(t) lc("$p2") mc("$p2") lwidth(thin))
            (connect `sharevar'_tsh_dfa_top1 year, ms(d) lc("$u4") mc("$u4") lpattern("-."))
                if inrange(year, 1960, 2019)
            ,
            $gpr ytitle("Top 1% `cmptname' Sh Total Wealth (%)") xtitle("")
            legend(order(2 "Equal" 1 "Baseline" 3 "SCF" 4 "DFA") 
                    region(lcolor(white) margin(tiny)) row(1))
            ysca(range(0 15)) ylab(0(2)14)
            xsize(5.5);
        graph_export `sharevar'_tshares_top1;

        twoway (connect `sharevar'_tsh_preferred_top01 year, ms(oh) lc("$u1") mc("$u1") lwidth(medthick))
            (connect `sharevar'_tsh_equ_top01 year,  ms(s) lc("$u3") mc("$u3")) 
            (connect `sharevar'_tsh_scf_top01 year,  ms(t) lc("$p2") mc("$p2") lwidth(thin))
                if inrange(year, 1960, 2019)
            ,
            $gpr ytitle("Top 0.1% `cmptname' Sh Total Wealth (%)") xtitle("")
            legend(order(2 "Equal" 1 "Baseline" 3 "SCF") 
                    region(lcolor(white) margin(tiny)) row(1))
            ysca(range(0 15)) ylab(0(2)14)
            xsize(5.5);
        graph_export `sharevar'_tshares_top01;

        twoway (connect `sharevar'_tsh_preferred_top001 year, ms(oh) lc("$u1") mc("$u1") lwidth(medthick))
            (connect `sharevar'_tsh_equ_top001 year,  ms(s) lc("$u3") mc("$u3")) 
            (connect `sharevar'_tsh_scf_top001 year, ms(t) lc("$p2") mc("$p2") lwidth(thin))
                if inrange(year, 1960, 2019)
            ,
            $gpr ytitle("Top 0.01% `cmptname' Sh Total Wealth (%)") xtitle("")
            legend(order(2 "Equal" 1 "Baseline" 3 "SCF") 
                    region(lcolor(white) margin(tiny)) row(1))
            ysca(range(0 15)) ylab(0(2)14)
            xsize(5.5);
        graph_export `sharevar'_tshares_top001;

        twoway (connect `sharevar'_tsh_preferred_p9099 year, ms(oh) lc("$u1") mc("$u1") lwidth(medthick))
            (connect `sharevar'_tsh_equ_p9099 year, ms(s) lc("$u3") mc("$u3"))
            (connect `sharevar'_tsh_scf_p9099 year, ms(t) lc("$p2") mc("$p2") lwidth(thin))
            (connect `sharevar'_tsh_dfa_p9099 year, ms(d) lc("$u4") mc("$u4") lpattern("-."))
                if inrange(year, 1960, 2019)
            ,
            $gpr ytitle("P90-99 `cmptname' Sh Total Wealth (%)") xtitle("")
            legend(order(2 "Equal" 1 "Baseline" 3 "SCF" 4 "DFA") 
                    region(lcolor(white) margin(tiny)) row(1))
            ysca(range(0 15)) ylab(0(2)14)
            xsize(5.5);
        graph_export `sharevar'_tshares_p9099;

        twoway (connect `sharevar'_tsh_preferred_bot90 year, ms(oh) lc("$u1") mc("$u1") lwidth(medthick))
            (connect `sharevar'_tsh_equ_bot90 year, ms(s) lc("$u3") mc("$u3"))
            (connect `sharevar'_tsh_scf_bot90 year, ms(t) lc("$p2") mc("$p2") lwidth(thin))
            (connect `sharevar'_tsh_dfa_bot90 year, ms(d) lc("$u4") mc("$u4") lpattern("-."))
                if inrange(year, 1960, 2019)
            ,
            $gpr ytitle("P0-90 `cmptname' Sh Total Wealth (%)") xtitle("")
            legend(order(2 "Equal" 1 "Baseline" 3 "SCF" 4 "DFA") 
                    region(lcolor(white) margin(tiny)) row(1))
            ysca(range(0 15)) ylab(0(2)14)
            xsize(5.5);
        graph_export `sharevar'_tshares_bot90;
        #delimit cr
    }

    foreach sharevar in hwfix ccorw pthru hwhou hwpen  {
    
        #delimit ;
        local cmptname = cond("`sharevar'" == "hwfix", "Fixed Income", 
                         cond("`sharevar'" == "hwequ", "Public Equity", 
                         cond("`sharevar'" == "hwbus", "Private Business", 
                         cond("`sharevar'" == "pthru", "Pass-Through Business", 
                         cond("`sharevar'" == "ccorw", "C-corporation Equity", 
                         cond("`sharevar'" == "hwpen", "Pensions", 
                         cond("`sharevar'" == "hwhou", "Housing", 
                                "Pensions, Housing, Other")))))));

        twoway (connect `sharevar'_preferred_top1 year, ms(oh) lc("$u1") mc("$u1") lwidth(medthick))
            (connect `sharevar'_equ_top1 year, ms(s) lc("$u3") mc("$u3")) 
            (connect `sharevar'_scf_top1 year, ms(t) lc("$p2") mc("$p2") lwidth(thin))
            (connect `sharevar'_dfa_top1 year, ms(d) lc("$u4") mc("$u4") lpattern("-."))
                if inrange(year, 1960, 2019)
            ,
            $gpr ytitle("Top 1% `cmptname' Levels (\$ T, 2016)") xtitle("")
            legend(order(2 "Equal" 1 "Baseline" 3 "SCF" 4 "DFA") 
                    region(lcolor(white) margin(tiny)) row(1))
            ysca(range(0 15)) ylab(0(2)14) 
            xsize(5.5);
        graph_export `sharevar'_total_top1;

        twoway (connect `sharevar'_preferred_top01 year, ms(oh) lc("$u1") mc("$u1") lwidth(medthick))
            (connect `sharevar'_equ_top01 year,  ms(s) lc("$u3") mc("$u3")) 
            (connect `sharevar'_scf_top01 year,  ms(t) lc("$p2") mc("$p2") lwidth(thin))
                if inrange(year, 1960, 2019)
            ,
            $gpr ytitle("Top 0.1% `cmptname' Levels (\$ T, 2016)") xtitle("")
            legend(order(2 "Equal" 1 "Baseline" 3 "SCF") 
                    region(lcolor(white) margin(tiny)) row(1))
            ysca(range(0 15)) ylab(0(2)14) 
            xsize(5.5);
        graph_export `sharevar'_total_top01;

        twoway (connect `sharevar'_preferred_top001 year, ms(oh) lc("$u1") mc("$u1") lwidth(medthick))
            (connect `sharevar'_equ_top001 year,  ms(s) lc("$u3") mc("$u3")) 
            (connect `sharevar'_scf_top001 year, ms(t) lc("$p2") mc("$p2") lwidth(thin))
                if inrange(year, 1960, 2019)
            ,
            $gpr ytitle("Top 0.01% `cmptname' Levels (\$ T, 2016)") xtitle("")
            legend(order(2 "Equal" 1 "Baseline" 3 "SCF") 
                    region(lcolor(white) margin(tiny)) row(1))
            ysca(range(0 15)) ylab(0(2)14) 
            xsize(5.5);
        graph_export `sharevar'_total_top001;

        twoway (connect `sharevar'_preferred_p9099 year, ms(oh) lc("$u1") mc("$u1") lwidth(medthick))
            (connect `sharevar'_equ_p9099 year, ms(s) lc("$u3") mc("$u3"))
            (connect `sharevar'_scf_p9099 year, ms(t) lc("$p2") mc("$p2") lwidth(thin))
            (connect `sharevar'_dfa_p9099 year, ms(d) lc("$u4") mc("$u4") lpattern("-."))
                if inrange(year, 1960, 2019)
            ,
            $gpr ytitle("P90-99 `cmptname' Levels (\$ T, 2016)") xtitle("")
            legend(order(2 "Equal" 1 "Baseline" 3 "SCF" 4 "DFA") 
                    region(lcolor(white) margin(tiny)) row(1))
            ysca(range(0 15)) ylab(0(2)14) 
            xsize(5.5);
        graph_export `sharevar'_total_p9099;

        twoway (connect `sharevar'_preferred_bot90 year, ms(oh) lc("$u1") mc("$u1") lwidth(medthick))
            (connect `sharevar'_equ_bot90 year, ms(s) lc("$u3") mc("$u3"))
            (connect `sharevar'_scf_bot90 year, ms(t) lc("$p2") mc("$p2") lwidth(thin))
            (connect `sharevar'_dfa_bot90 year, ms(d) lc("$u4") mc("$u4") lpattern("-."))
                if inrange(year, 1960, 2019)
            ,
            $gpr ytitle("P0-90 `cmptname' Levels (\$ T, 2016)") xtitle("")
            legend(order(2 "Equal" 1 "Baseline" 3 "SCF" 4 "DFA") 
                    region(lcolor(white) margin(tiny)) row(1))
            ysca(range(0 15)) ylab(0(2)14) 
            xsize(5.5);
        graph_export `sharevar'_total_bot90;
        #delimit cr
    }/*%>*/

end/*%>*/

capture program drop graph_uncertainty/*%<*/
program define graph_uncertainty

    tempfile fig1
    insheet using topshares_adjusted.csv, comma clear
    keep year top*_equrtrns_v3 top*_preferred
    save `fig1'

    tempfile forbes
    load_analysis_data scfharmony
    gen forbes_share = forbeswlth / total_wealth
    keep year forbes_share
    save `forbes'

    * Produced by appendix program scf_stderrors_build()
    tempfile scf_ses
    use $dumpdir/scftopshares20220112.dta, clear

    ds t*share_*
    local checklist = "`r(varlist)'"

    rename t*share_* t*share_*_check

    merge 1:1 year using $dumpdir/scftopshares_stderrors20220112.dta, assert(3) nogen
    merge 1:1 year using `forbes', keep(3) nogen

    foreach v of varlist t*share* {
        replace `v' = `v' + forbes_share
        replace `v' = 100 * `v'
    }

    save `scf_ses'

    foreach series in "wealth_share_top01" "wealth_share_top001" "wealth_share_top1" {
        use "$inputs/uncertainty_numbers", clear

        *keep the series we are plotting
        keep `series' year 
        replace `series' = 100 * `series'

        *remove the top and bottom 2.5% to estimate 95% interval
        sort year
        by year: egen p2pt5 = pctile(`series'), p(2.5)
        by year: egen p97pt5 = pctile(`series'), p(97.5)
        by year: egen mean = mean(`series')

        drop if `series' > p97pt5 | `series' < p2pt5

        merge m:1 year using `scf_ses', keep(1 2 3) nogen
        merge m:1 year using `fig1', keep(3) nogen

        /************************************************************************START*
        Create plot
        *************************************************************************END*/

        local series_name = cond("`series'" == "wealth_share_top01", "Top 0.1% Wealth Share", ///
                            cond("`series'" == "taxbond_share_top01", "Top 0.1% Taxable Fixed Income Wealth as" "Share of Net Household Wealth (%)", ///
                            cond("`series'" == "ccorp_share_top01", "C Corp share", ///
                                "ERROR")))

        local yscale =     cond("`series'" == "wealth_share_top01", "ysca(range(5 20)) ylab(6(2)20)", ///
                            cond("`series'" == "wealth_share_top001", "ysca(range(1 11)) ylab(2(2)10)", ///
                            cond("`series'" == "wealth_share_top1", "ysca(range(18 41)) ylab(20(5)40)", ///
                            cond("`series'" == "taxbond_share_top01", "", ///
                            cond("`series'" == "ccorp_share_top01", "", ///
                                "ERROR")))))

        local scfvar = cond("`series'" == "wealth_share_top01", "t01share_pref", ///
                            cond("`series'" == "wealth_share_top001", "t001share_pref", ///
                            cond("`series'" == "wealth_share_top1", "t1share_pref", "")))

        local capvar = cond("`series'" == "wealth_share_top01", "top01", ///
                            cond("`series'" == "wealth_share_top001", "top001", ///
                            cond("`series'" == "wealth_share_top1", "top1", "")))

        *Create new counter variable after removing the top and bottom 2.5%
        bysort year: gen counter = _n
        
        *reset local
        local graphs
        forvalues  i = 1(1)100 {
            local graphs `graphs' (line `series' year if counter == `i', color(gray*.7) lwidth(vvvthin))
        }
        gsort counter year
        #delimit ;
        graph twoway `graphs' 
            (connect mean year if counter == 1, color("$u1") ms(none))
            (connect `scfvar' year if counter == 1, ms(t) color("$p2"))
            (connect `scfvar'_ci95_lb year if counter == 1, ms(none) color("$p2") lp("_"))
            (connect `scfvar'_ci95_ub year if counter == 1, ms(none) color("$p2") lp("_"))
            (connect `capvar'_pref year if counter == 1 & year >= 2001, ms(o) color("$u1") lw(medthick)) 
            (connect `capvar'_equrtrns_v3 year if counter == 1, ms(s) color("$u3")) 
            , 
            legend(order(106 "Equal Returns"
                         105 "Baseline" 101 "CMD 3-Tier" 100 "3-Tier Bstrap"
                         102 "Harm. SCF" 103 "SCF 95% CI") symxsize(*.4) row(2) region(lc(white)))
            ytitle("Share of Total Household Wealth (%)") xtitle("") xsize(4.5)
            `yscale'
            xlab(1960(20)2020) 
         $gpr;
         #delimit cr

        graph_export "`series'_imputed_interval"
    }

end/*%>*/

capture program drop graph_imputation/*%<*/
program define graph_imputation

    tempfile sz2020
    insheet using sz2020_tu_top01.csv, comma clear
    replace top01_sz20 = 100 * top01_sz20
    save `sz2020'

    tempfile szztu
    load_analysis_data szz_tu
    keep if inlist(group, "All", "P99.9-100")
    sort year group
    by year: assert _N == 2
    replace group = cond(group == "All", "_all", "_top01")
    keep year group hweal_preferred
    reshape wide hweal_preferred, i(year) j(group, string)
    gen top01_pref = (hweal_preferred_top01 / hweal_preferred_all) * 100
    save `szztu'

    tempfile fig1
    insheet using topshares_adjusted.csv, comma clear
    keep year top*_equrtrns_v3 top*_preferred
    save `fig1'

    tempfile fig1a
    insheet using top01shares_unadjusted.csv, comma clear
    keep year top01_sz
    save `fig1a'

    tempfile forbes
    load_analysis_data scfharmony_v3
    gen forbes_share = forbeswlth / total_wealth
    keep year forbes_share
    save `forbes'

    * Produced by appendix program scf_stderrors_build()
    tempfile scf_ses
    use $dumpdir/scftopshares20210720.dta, clear

    ds t*share_*
    local checklist = "`r(varlist)'"

    rename t*share_* t*share_*_check

    merge 1:1 year using $dumpdir/scftopshares_stderrors20210720.dta, assert(3) nogen
    merge 1:1 year using `forbes', keep(3) nogen

    foreach v of varlist t*share* {
    replace `v' = `v' + forbes_share
    replace `v' = 100 * `v'
    }

    save `scf_ses'

    tempfile fig2
    insheet using aggwealth.csv, comma clear
    gen hweal_scale = (hweal_preferred - 1e6 * ttpen + ttpen_preferred) / hweal_preferred
    gen pen_scale = ttpen_preferred / (1e6 * ttpen)
    keep year hweal_scale pen_scale
    save `fig2'

    tempfile fig8
    insheet using social_security_top01.csv, comma clear
    keep year top01sh_ssw_shv
    save `fig8'

    * Created by graph_comparison_trends
    use imputation_preferred.dta, clear
    merge 1:1 year using `fig1', nogen
    merge 1:1 year using `fig2', nogen
    merge 1:1 year using `fig8', nogen
    merge 1:1 year using `scf_ses', keep(1 2 3) nogen

    * 1. szz (ie pref)
    #delimit ;
    gen spec_pref_chk = top01_taxbond_pref_sh_aggwlth + top01_divw_pref_sh_aggwlth +
                    top01_pthru_base_sh_aggwlth + top01_hwhou_pref_sh_aggwlth + 
                    top01_hwpen_base_sh_aggwlth + top01_oth_base_sh_aggwlth;
    #delimit cr
    rename top01_preferred spec_pref
    * 1b szz cmd (ie cmd 3 tier)
    gen spec_1b = spec_pref - top01_taxbond_pref_sh_aggwlth + top01_taxbond_cmd3_sh_aggwlth
    * [1d] no fa agg
    * 2. szz pref
    * 2b. .75 divs
    gen spec_2b = spec_pref - top01_divw_pref_sh_aggwlth + top01_divw_25kg_sh_aggwlth
    * 3. szz pref
    * 3b unscaled pass-through
    gen spec_3b = spec_pref - top01_pthru_base_sh_aggwlth + top01_pthru_pref_sh_aggwlth
    * 3c no labor adj
    gen spec_3c = spec_pref - top01_pthru_base_sh_aggwlth + top01_pthru_full_sh_aggwlth
    * 3d no liquid
    gen spec_3d = spec_3b + (.1/.9) * top01_pthru_pref_sh_aggwlth 
    * 3f equal including missing
    gen spec_3f = spec_pref - top01_pthru_base_sh_aggwlth + top01_pthru_equM_sh_aggwlth 
    * 4. pref
    * [4b] with unfunded db
    gen spec_4b = (spec_pref - top01_hwpen_base_sh_aggwlth) / hweal_scale
    replace spec_4b = spec_4b + top01_hwpen_pref_sh_aggwlth 
    * [4c] with soc security
    gen spec_4c = top01sh_ssw_shv
    * 5 pref
    * 5b equal returns
    gen spec_5b = spec_pref - top01_hwhou_pref_sh_aggwlth + top01_hwhou_equ_sh_aggwlth

    keep if year >= 1966 & year <= 2019
    gsort year
    #delimit ;
    graph twoway 
    (connect spec_1b year, msym(oh) msize(small) color(gray*.7) )
    (connect spec_2b year, msym(oh) msize(small) color(gray*.7) )
    (connect spec_3b year, msym(oh) msize(small) color(gray*.7) )
    (connect spec_3c year, msym(oh) msize(small) color(gray*.7) )
    (connect spec_3d year, msym(oh) msize(small) color(gray*.7) )
    (connect spec_3f year, msym(oh) msize(small) color(gray*.7) )
    (connect spec_4b year, msym(oh) msize(small) color(gray*.7) )
    (connect spec_4c year, msym(s) msize(small) color("$f1") )
    (connect spec_5b year, msym(oh) msize(small) color(gray*.7) )
    (connect spec_pref year, ms(oh) msize(small) color("$u1") lw(medthin)) 
    (connect t01share_pref year, ms(t) color("$p2"))
    (connect t01share_pref_ci95_lb year, ms(none) color("$p2") lp("_"))
    (connect t01share_pref_ci95_ub year, ms(none) color("$p2") lp("_"))
    (connect top01_equrtrns_v3 year, ms(s) color("$u3")) 
    , 
    legend(order(14 "Equal Returns"
                 10 "Baseline" 2 "Perturbation"
                 11 "Harm. SCF" 12 "SCF 95% CI" 8 "Base. w/Soc Sec") symxsize(*.4) row(1) region(lc(white)))
    ytitle("Share of Total Household Wealth (%)") xtitle("") xsize(8)
    ysca(range(5 20)) ylab(5(2.5)20)
    xlab(1965(5)2015) 
    $gpr;
    #delimit cr

    graph_export "imputation_interval"

end/*%>*/

capture program drop table_portfolioshares/*%<*/
program define table_portfolioshares
/*******************************************************************************
    Portfolio Shares by Wealth Group. Shows portfolio shares under 
        preferred portfolio construction and rankings.
*******************************************************************************/

    /***************************************************************************
        Prep to build head of the table (e.g. title, etc.) via listtex, 
            which doesn't depend on data
    ***************************************************************************/

    clear
    qui set obs 1

    qui gen tab = "\begin{tabular}{lcccccc}" in 1

    qui gen hline = "\hline" in 1
    qui gen cline = "\cmidrule(lr){2-7}" in 1
    gen midrule = "\midrule" in 1

    qui gen c = " " in 1

    qui gen title1 = " & \multicolumn{6}{c}{Portfolio shares (Baseline)}" in 1

    #delimit ;
    qui gen title2 = "Wealth group & Fixed Income & C-corporation Equity & 
                      Pass-through Business & Housing & Pensions & Other " in 1;
    #delimit cr

    local tabname = "portfolioshares"

    *Output
    qui listtex tab if _n == 1 using "`tabname'.tex", replace rstyle(none)
    qui listtex midrule if _n == 1, appendto("`tabname'.tex") rstyle(none)
    qui listtex title1 if _n == 1, appendto("`tabname'.tex") rstyle(tabular)
    qui listtex cline if _n == 1, appendto("`tabname'.tex") rstyle(none)
    qui listtex title2 if _n == 1, appendto("`tabname'.tex") rstyle(tabular)
    qui listtex midrule if _n==1 , appendto("`tabname'.tex") rstyle(none)    

    /***************************************************************************
        Cycling through ranking specifications, make data-dependent 
            parts of table, namely wealth totals, levels, thresholds, and group
            counts under each ranking specification: SZ (baseline) and preferred
    ***************************************************************************/
    
        /***********************************************************************
            Load in tax data using dina_collapse program for year 2016;
                build portfolio categories according to proper capitalization 
                specification
        ***********************************************************************/

    load_analysis_data szz
    keep if year == 2016

    keep w${preferred_defn_late}_group group hweal$preferred_defn_late *_pref *_base

    egen hweal_check = rowtotal(hwfix_base pthru_base hwpen_base hwhou_base ccorw_base hwoth_base)
    assert inrange(hweal_check / hweal$preferred_defn_late, 0.99, 1.01)
    drop hweal_check

    rename *_base *

        /***********************************************************************
            Drop unnecessary wealth groups rows; keep only mutually 
                inclusive rows. Ensure 5 wealth groups after dropping these 
                wealth groups. 
        ***********************************************************************/

    sort w${preferred_defn_late}_group

    qui keep if inrange(w${preferred_defn_late}_group, 1, 11)

    assert _N == 11

        /***********************************************************************
            Compute wealth shares and average wealth in groups
        ***********************************************************************/

    foreach portfvar of varlist pthru ccorw hwfix hwhou hwpen hwoth {
        qui gen sh_`portfvar'  = 100 * (`portfvar' / hweal$preferred_defn_late)
        drop `portfvar'
    }
    drop hweal$preferred_defn_late

        /***********************************************************************
            Reformat variables in preparation to put in table
        ***********************************************************************/

    * Label wealth groups
    gen grplab = cond(group == "All", "Full population", ///
                    cond(group == "P90-100", "Top 10\%", ///
                    cond(group == "P99-100", "Top 1\%", ///
                    cond(group == "P99.9-100", "Top 0.1\%", ///
                    cond(group == "P99.99-100", "Top 0.01\%", ///
                    cond(group == "P99.999-100", "Top 0.001\%", ///
                    cond(group == "P0-90", "Bottom 90\%", ///
                    cond(group == "P90-99", "Top 10-1\%", ///
                    cond(group == "P99-99.9", "Top 1-0.1\%", ///
                    cond(group == "P99.9-99.99", "Top 0.1-0.01\%", ///
                    "Top 0.01-0.001\%"))))))))))

    export delim using portfolioshares.csv, replace

    * Converting to strings
    qui format sh* %12.1f
    qui tostring sh*, replace usedisplayformat force

    foreach sharevar of varlist sh* {
        qui replace `sharevar' = `sharevar' + "\%"
    }

        /***********************************************************************
            Write data-dependent part of table under ranking in question
                via listtex
        ***********************************************************************/

    qui gen panelA = "\multicolumn{7}{l}{\textbf{Panel A. Top wealth groups}}" in 1
    qui gen panelB = "\multicolumn{7}{l}{\textbf{Panel B. Intermediate wealth groups}}" in 1
    
    qui gen hline = "\hline" in 1
    qui gen c = " " in 1
    qui gen end = "\end{tabular}" in 1
    gen vsp = "\vspace{-.1in}" in 1

    qui listtex panelA if _n == 1 , appendto("`tabname'.tex") rstyle(tabular)
    
    qui listtex grplab sh_hwfix sh_ccorw sh_pthru sh_hwhou sh_hwpen sh_hwoth ///
        if w${preferred_defn_late}_group <= 6, appendto("`tabname'.tex") ///
        rstyle(tabular)
    
    listtex vsp if _n == 1, appendto("`tabname'.tex") rstyle(tabular)   
    qui listtex panelB if _n == 1 , appendto("`tabname'.tex") rstyle(tabular)
    
    qui listtex grplab sh_hwfix sh_ccorw sh_pthru sh_hwhou sh_hwpen sh_hwoth ///
        if w${preferred_defn_late}_group > 6, appendto("`tabname'.tex") ///
        rstyle(tabular)

    /***************************************************************************
        Close out table
    ***************************************************************************/

    qui listtex hline if _n == 1, appendto("`tabname'.tex") rstyle(none)  
    qui listtex end if _n == 1, appendto("`tabname'.tex") rstyle(none)

end/*%>*/

capture program drop setup_suppA/*%<*/
program define setup_suppA

    tempfile outdata
    
    ****************************************************************************
    * Panel A: Top Wealth Shares and Per-Capita Wealth in 2016
    ****************************************************************************
    ****************************************************************************
    * 1. Baseline
    ****************************************************************************
    load_taxdata, rankspec($preferred_defn_late) year(2016)

    #delimit ;
    keep if inlist(group, "All", "P99-100", "P99.9-100", "P99.99-100", "P99.999-100");

    replace group = cond(group == "P99.999-100", "_top0001",
                    cond(group == "P99.99-100", "_top001", 
                    cond(group == "P99.9-100", "_top01", 
                    cond(group == "P99-100", "_top1", 
                    cond(group == "All", "_total", "")))));
    drop if group == "";
    #delimit cr
    
    keep year group n hweal$preferred_defn_late
    rename hweal$preferred_defn_late hweal
    reshape wide hweal n, i(year) j(group, string)

    foreach topgrp in top1 top01 top001 top0001 {
        gen `topgrp'_share = (hweal_`topgrp' / hweal_total) * 100
        gen `topgrp'_percap = hweal_`topgrp' / n_`topgrp'
    }
    gen spec = "baseline"
    gen speclab = "1. Baseline"

    sum hweal_total
    local aggbase = `r(mean)'
    save `outdata'

    tempfile ngrp
    keep year n_*
    save `ngrp'
    ****************************************************************************
    * 2. Equal returns
    ****************************************************************************
    load_taxdata, rankspec($sz_wlth_defn_late_v3) year(2016)

    #delimit ;
    keep if inlist(group, "All", "P99-100", "P99.9-100", "P99.99-100", "P99.999-100");

    replace group = cond(group == "P99.999-100", "_top0001",
                    cond(group == "P99.99-100", "_top001", 
                    cond(group == "P99.9-100", "_top01", 
                    cond(group == "P99-100", "_top1", 
                    cond(group == "All", "_total", "")))));
    drop if group == "";
    #delimit cr
    
    keep year group n hweal$sz_wlth_defn_late_v3
    rename hweal$sz_wlth_defn_late_v3 hweal
    reshape wide hweal n, i(year) j(group, string)

    foreach topgrp in top1 top01 top001 top0001 {
        gen `topgrp'_share = (hweal_`topgrp' / hweal_total) * 100
        gen `topgrp'_percap = hweal_`topgrp' / n_`topgrp'
    }
    gen spec = "equalrets"
    gen speclab = "2. Equal Returns"

    append using `outdata'
    save `outdata', replace
    ****************************************************************************
    * 3. SZ 2020
    ****************************************************************************
    load_analysis_data sz2020_top
    keep if year == 2016
    rename top*_pszes20 top*_share
    drop top10_share

    * By construction, should match revised SZ
    foreach topgrp in top1 top01 top001 top0001 {
        gen hweal_`topgrp' = .01 * `topgrp'_share * `aggbase' 
    }
    gen hweal_total = `aggbase'

    gen spec = "revisedsz"
    gen speclab = "3. Revised Saez Zucman (2020)"

    append using `outdata'
    foreach topgrp in top1 top01 top001 top0001 total {
        egen n_`topgrp'HB = mean(n_`topgrp') if spec == "baseline"
        egen n_`topgrp'H = max(n_`topgrp'HB)
        replace n_`topgrp' = n_`topgrp'H if spec == "revisedsz"
        drop n_`topgrp'H*
    }

    foreach topgrp in top1 top01 top001 top0001 {
        replace `topgrp'_percap = hweal_`topgrp' / n_`topgrp' if spec == "revisedsz"
    }

    save `outdata', replace
    ****************************************************************************
    * 4. PSZ 2018
    ****************************************************************************
    load_taxdata, rankspec($sz_wlth_defn) year(2016)

    #delimit ;
    keep if inlist(group, "All", "P99-100", "P99.9-100", "P99.99-100", "P99.999-100");

    replace group = cond(group == "P99.999-100", "_top0001",
                    cond(group == "P99.99-100", "_top001", 
                    cond(group == "P99.9-100", "_top01", 
                    cond(group == "P99-100", "_top1", 
                    cond(group == "All", "_total", "")))));
    drop if group == "";
    #delimit cr
    
    keep year group n hweal$sz_wlth_defn
    rename hweal$sz_wlth_defn hweal
    reshape wide hweal n, i(year) j(group, string)

    foreach topgrp in top1 top01 top001 top0001 {
        gen `topgrp'_share = (hweal_`topgrp' / hweal_total) * 100
        gen `topgrp'_percap = hweal_`topgrp' / n_`topgrp'
    }
    gen spec = "psz"
    gen speclab = "4. Piketty Saez Zucman (2018)"

    append using `outdata'
    save `outdata', replace
    ****************************************************************************
    * 5, 6, 7, 8, 9, 14: Alternative Aggregates
    ****************************************************************************
    * 5: 1 with unscaled pass-through
    * 6: 5 with missing pass-through
    * 7: 6 with unfunded pensions
    * 8: 7 with debt adjustments
    * 9: 8 with no student debt
    * 14: 2 with bond split (robustness)
    * 15: 1 with trimmed interest rates
    * 16: 8 with trimmed interest rates
    ****************************************************************************
    foreach rankspec in 3021 3022 3023 3024 3025 3018 3026 3027 {
        load_taxdata, rankspec(`rankspec') year(2016)

        #delimit ;
        keep if inlist(group, "All", "P99-100", "P99.9-100", "P99.99-100", "P99.999-100");

        replace group = cond(group == "P99.999-100", "_top0001",
                        cond(group == "P99.99-100", "_top001", 
                        cond(group == "P99.9-100", "_top01", 
                        cond(group == "P99-100", "_top1", 
                        cond(group == "All", "_total", "")))));
        drop if group == "";
        #delimit cr
        
        keep year group n hweal`rankspec'
        rename hweal`rankspec' hweal
        reshape wide hweal n, i(year) j(group, string)

        foreach topgrp in top1 top01 top001 top0001 {
            gen `topgrp'_share = (hweal_`topgrp' / hweal_total) * 100
            gen `topgrp'_percap = hweal_`topgrp' / n_`topgrp'
        }
        gen spec = "alt`rankspec'"

        if (`rankspec' == 3021) {
            gen speclab = "5. Spec \#1 with Unscaled Pass-through"
        }
        else if (`rankspec' == 3022) {
            gen speclab = "6. Spec \#5 with Missing Pass-through"
        }
        else if (`rankspec' == 3023) {
            gen speclab = "7. Spec \#6 with Unfunded Pensions"
        }
        else if (`rankspec' == 3024) {
            gen speclab = "8. Spec \#7 with Debt Adjustments"
        }
        else if (`rankspec' == 3025) {
            gen speclab = "9. Spec \#8 with No Student Debt"
        }
        else if (`rankspec' == 3018) {
            gen speclab = "14. Spec \#2 with Bond Split"
        }
        else if (`rankspec' == 3026) {
            gen speclab = "15. Spec \#1 with Low Boutique"
        }
        else if (`rankspec' == 3027) {
            gen speclab = "16. Spec \#8 with Low Boutique"
        }
        append using `outdata'
        save `outdata', replace
    }
    ****************************************************************************
    * 10, 11, 12, 13, 17: Forbes specs
    ****************************************************************************
    * 10: Baseline with BHV
    * 11: Baseline with Replace
    * 12: 6 with BHV
    * 13: 8 with BHV
    * 17: 13 with BHV
    ****************************************************************************
    foreach rankspec in 3016 3022 3024 3027 {
        load_analysis_data szz_forbes `rankspec'
        keep if year == 2016
        merge 1:1 year using `ngrp', keep(3) nogen

       rename hweal*_adjusted hweal*
        foreach topgrp in top1 top01 top001 top0001 {
            gen `topgrp'_share = `topgrp'share_adjusted 
            gen `topgrp'_percap = hweal_`topgrp' / n_`topgrp'
        }
        gen spec = "forbes`rankspec'"

        keep year spec hweal_total *_share *_percap n_total
        if ("`rankspec'" == "3016") {
            gen speclab = "10. Spec \#1 with BHV Reweighting"
        }
        else if ("`rankspec'" == "3022") {
            gen speclab = "12. Spec \#6 with BHV Reweighting"
        }
        else if ("`rankspec'" == "3024") {
            gen speclab = "13. Spec \#8 with BHV Reweighting"
        }
        else if ("`rankspec'" == "3027") {
            gen speclab = "17. Spec \#13 with Low Boutique"
        }
        append using `outdata'
        save `outdata', replace
    }

    load_analysis_data szz_forbes 3016
    keep if year == 2016
    merge 1:1 year using `ngrp', keep(3) nogen

    rename hweal*_replaced hweal*
    foreach topgrp in top1 top01 top001 top0001 {
        gen `topgrp'_share = `topgrp'share_replaced 
        gen `topgrp'_percap = hweal_`topgrp' / n_`topgrp'
    }
    gen spec = "forbes3016replace"

    keep year spec hweal_total *_share *_percap n_total
    gen speclab = "11. Spec \#1 with Replace Top 400"
    append using `outdata'
    save `outdata', replace

    use `outdata', clear

end/*%>*/

capture program drop setup_suppB/*%<*/
program define setup_suppB

    tempfile outdata
    
    ****************************************************************************
    * Panel B: Top Portfolio Shares in 2016 (%)
    ****************************************************************************
    ****************************************************************************
    * 1. Baseline
    ****************************************************************************
    load_taxdata, rankspec($preferred_defn_late) year(2016)

    #delimit ;
    keep if inlist(group, "All", "P99-100", "P99.9-100", "P99.99-100", "P99.999-100");

    replace group = cond(group == "P99.999-100", "_top0001",
                    cond(group == "P99.99-100", "_top001", 
                    cond(group == "P99.9-100", "_top01", 
                    cond(group == "P99-100", "_top1", 
                    cond(group == "All", "_total", "")))));
    drop if group == "";
    #delimit cr
    
    keep year group hweal$preferred_defn_late hwfix_base pthru_base ccorw_base 
    rename hweal$preferred_defn_late hweal
    rename *_base *
    reshape wide hweal hwfix pthru ccorw, i(year) j(group, string)

    foreach topgrp in top1 top01 top001 top0001 {
        gen `topgrp'_sh_hwfix = (hwfix_`topgrp' / hweal_`topgrp') * 100
        gen `topgrp'_sh_pthru = (pthru_`topgrp' / hweal_`topgrp') * 100
        gen `topgrp'_sh_ccorw = (ccorw_`topgrp' / hweal_`topgrp') * 100
    }
    gen spec = "baseline"
    gen speclab = "1. Baseline"

    sum hweal_total
    local aggbase = `r(mean)'
    save `outdata'
    ****************************************************************************
    * 2. Equal returns
    ****************************************************************************
    load_taxdata, rankspec($sz_wlth_defn_late_v3) year(2016)

    #delimit ;
    keep if inlist(group, "All", "P99-100", "P99.9-100", "P99.99-100", "P99.999-100");

    replace group = cond(group == "P99.999-100", "_top0001",
                    cond(group == "P99.99-100", "_top001", 
                    cond(group == "P99.9-100", "_top01", 
                    cond(group == "P99-100", "_top1", 
                    cond(group == "All", "_total", "")))));
    drop if group == "";
    #delimit cr
    
    keep year group hweal$sz_wlth_defn_late_v3 hwfix_equ pthru_equ ccorw_equ
    rename hweal$sz_wlth_defn_late_v3 hweal
    rename *_equ *
    reshape wide hweal hwfix pthru ccorw, i(year) j(group, string)

    foreach topgrp in top1 top01 top001 top0001 {
        gen `topgrp'_sh_hwfix = (hwfix_`topgrp' / hweal_`topgrp') * 100
        gen `topgrp'_sh_pthru = (pthru_`topgrp' / hweal_`topgrp') * 100
        gen `topgrp'_sh_ccorw = (ccorw_`topgrp' / hweal_`topgrp') * 100
    }
    gen spec = "equalrets"
    gen speclab = "2. Equal Returns"

    append using `outdata'
    save `outdata', replace
    ****************************************************************************
    * 3. SZ 2020
    ****************************************************************************
    load_analysis_data sz2020_port
    keep if year == 2016

    foreach topgrp in top1 top01 top001 {
        gen `topgrp'_sh_hwfix = (fix_`topgrp' / wlth_`topgrp') * 100
        gen `topgrp'_sh_pthru = .
        gen `topgrp'_sh_ccorw = .
    }

    gen spec = "revisedsz"
    gen speclab = "3. Revised Saez Zucman (2020)"

    append using `outdata'
    save `outdata', replace
    ****************************************************************************
    * 4. PSZ 2018
    ****************************************************************************
    load_taxdata, rankspec($sz_wlth_defn) year(2016)

    #delimit ;
    keep if inlist(group, "All", "P99-100", "P99.9-100", "P99.99-100", "P99.999-100");

    replace group = cond(group == "P99.999-100", "_top0001",
                    cond(group == "P99.99-100", "_top001", 
                    cond(group == "P99.9-100", "_top01", 
                    cond(group == "P99-100", "_top1", 
                    cond(group == "All", "_total", "")))));
    drop if group == "";
    #delimit cr
    
    keep year group hweal$sz_wlth_defn hwfix_equ pthru_equ ccorw_equ
    rename hweal$sz_wlth_defn hweal
    rename *_equ *
    reshape wide hweal hwfix pthru ccorw, i(year) j(group, string)

    foreach topgrp in top1 top01 top001 top0001 {
        gen `topgrp'_sh_hwfix = (hwfix_`topgrp' / hweal_`topgrp') * 100
        gen `topgrp'_sh_pthru = (pthru_`topgrp' / hweal_`topgrp') * 100
        gen `topgrp'_sh_ccorw = (ccorw_`topgrp' / hweal_`topgrp') * 100
    }
    gen spec = "psz"
    gen speclab = "4. Piketty Saez Zucman (2018)"

    append using `outdata'
    save `outdata', replace
    ****************************************************************************
    * 5, 6, 7, 8, 9: Alternative Aggregates
    ****************************************************************************
    * 5: 1 with unscaled pass-through
    * 6: 5 with missing pass-through
    * 7: 6 with unfunded pensions
    * 8: 7 with debt adjustments
    * 9: 8 with no student debt
    * 14: 2 with bond split (robustness)
    * 15: 1 with trimmed interest rates
    * 16: 8 with trimmed interest rates
    ****************************************************************************
    foreach rankspec in 3021 3022 3023 3024 3025 3018 3026 3027 {
        load_taxdata, rankspec(`rankspec') year(2016)

        #delimit ;
        keep if inlist(group, "All", "P99-100", "P99.9-100", "P99.99-100", "P99.999-100");

        replace group = cond(group == "P99.999-100", "_top0001",
                        cond(group == "P99.99-100", "_top001", 
                        cond(group == "P99.9-100", "_top01", 
                        cond(group == "P99-100", "_top1", 
                        cond(group == "All", "_total", "")))));
        drop if group == "";
        #delimit cr
        
        * Note hwfix_base = hwfix_pref and ccorw_base = ccorw_pref
        if (`rankspec' == 3021) {
            keep year group hweal`rankspec' hwfix_pref pthru_pref ccorw_pref ///
                missing_scorp missing_pship
            replace pthru_pref = pthru_pref - missing_scorp - missing_pship
            rename *_pref *
            drop missing_scorp missing_pship
        }
        else if inlist(`rankspec', 3022, 3023, 3024, 3025) {
            keep year group hweal`rankspec' hwfix_pref pthru_pref ccorw_pref 
            rename *_pref *
        }
		else if inlist(`rankspec', 3026) {
            keep year group hweal`rankspec' hwfix_pref pthru_base ccorw_pref ///
                taxbond_info_trim taxbond_info
            replace hwfix_pref = hwfix_pref - taxbond_info + taxbond_info_trim
            rename *_pref *
            drop taxbond*
			rename *_base *
        }
        else if inlist(`rankspec', 3027) {
            keep year group hweal`rankspec' hwfix_pref pthru_pref ccorw_pref ///
                taxbond_info_trim taxbond_info
            replace hwfix_pref = hwfix_pref - taxbond_info + taxbond_info_trim
            rename *_pref *
            drop taxbond*
        }
        else if inlist(`rankspec', 3018) {
            keep year group hweal`rankspec' hwfix_equal pthru_equ ccorw_equal 
            rename *_equal *
            rename *_equ *
        }
        rename hweal`rankspec' hweal
        reshape wide hweal hwfix pthru ccorw, i(year) j(group, string)

        foreach topgrp in top1 top01 top001 top0001 {
            gen `topgrp'_sh_hwfix = (hwfix_`topgrp' / hweal_`topgrp') * 100
            gen `topgrp'_sh_pthru = (pthru_`topgrp' / hweal_`topgrp') * 100
            gen `topgrp'_sh_ccorw = (ccorw_`topgrp' / hweal_`topgrp') * 100
        }
        gen spec = "alt`rankspec'"

        if (`rankspec' == 3021) {
            gen speclab = "5. Spec \#1 with Unscaled Pass-through"
        }
        else if (`rankspec' == 3022) {
            gen speclab = "6. Spec \#5 with Missing Pass-through"
        }
        else if (`rankspec' == 3023) {
            gen speclab = "7. Spec \#6 with Unfunded Pensions"
        }
        else if (`rankspec' == 3024) {
            gen speclab = "8. Spec \#7 with Debt Adjustments"
        }
        else if (`rankspec' == 3025) {
            gen speclab = "9. Spec \#8 with No Student Debt"
        }
        else if (`rankspec' == 3018) {
            gen speclab = "14. Spec \#2 with Bond Split"
        }
        else if (`rankspec' == 3026) {
            gen speclab = "15. Spec \#1 with Low Boutique"
        }
        else if (`rankspec' == 3027) {
            gen speclab = "16. Spec \#8 with Low Boutique"
        }
        append using `outdata'
        save `outdata', replace
    }
    ****************************************************************************
    * 10, 11, 12, 13: Forbes specs
    ****************************************************************************
    * 10: Baseline with BHV
    * 11: Baseline with Replace
    * 12: 6 with BHV
    * 13: 8 with BHV
    * 17: 13 with BHV
    ****************************************************************************
    foreach rankspec in 3016 3022 3024 3027 {
        load_analysis_data szz_forbes `rankspec'
        keep if year == 2016

        if inlist(`rankspec', 3022, 3024, 3027) {
            keep year hweal_* hwfix_pref_*_adjusted pthru_pref_*_adjusted ccorw_pref_*_adjusted ///
                forbes400_hwfix_adjusted forbes400_pthru_adjusted forbes400_ccorw_adjusted
            rename *_adjusted *
            rename hwfix_pref_* hwfix_*
            rename pthru_pref_* pthru_*
            rename ccorw_pref_* ccorw_*
        }
        else if inlist(`rankspec', 3016) {
            keep year hweal_* hwfix_base_*_adjusted pthru_base_*_adjusted ccorw_base_*_adjusted ///
                forbes400_hwfix_adjusted forbes400_pthru_adjusted forbes400_ccorw_adjusted
            rename *_adjusted *
            rename hwfix_base_* hwfix_*
            rename pthru_base_* pthru_*
            rename ccorw_base_* ccorw_*
        }

        foreach topgrp in top1 top01 top001 top0001 {
            gen `topgrp'_sh_hwfix = ((hwfix_`topgrp' + forbes400_hwfix) / hweal_`topgrp') * 100
            gen `topgrp'_sh_pthru = ((pthru_`topgrp' + forbes400_pthru) / hweal_`topgrp') * 100
            gen `topgrp'_sh_ccorw = ((ccorw_`topgrp' + forbes400_ccorw) / hweal_`topgrp') * 100
        }

        gen spec = "forbes`rankspec'"

        keep year spec *_sh_* 
        if ("`rankspec'" == "3016") {
            gen speclab = "10. Spec \#1 with Reweighting"
        }
        else if ("`rankspec'" == "3022") {
            gen speclab = "12. Spec \#6 with Reweighting"
        }
        else if ("`rankspec'" == "3024") {
            gen speclab = "13. Spec \#8 with Reweighting"
        }
        else if ("`rankspec'" == "3027") {
            gen speclab = "17. Spec \#13 with Low Boutique"
        }
        append using `outdata'
        save `outdata', replace
    }

    load_analysis_data szz_forbes 3016
    keep if year == 2016

    keep year hweal_* hwfix_base_*_replaced pthru_base_*_replaced ccorw_base_*_replaced ///
        forbes400_hwfix_replaced forbes400_pthru_replaced forbes400_ccorw_replaced
    rename *_replaced *
    rename hwfix_base_* hwfix_*
    rename pthru_base_* pthru_*
    rename ccorw_base_* ccorw_*

    foreach topgrp in top1 top01 top001 top0001 {
        gen `topgrp'_sh_hwfix = ((hwfix_`topgrp' + forbes400_hwfix) / hweal_`topgrp') * 100
        gen `topgrp'_sh_pthru = ((pthru_`topgrp' + forbes400_pthru) / hweal_`topgrp') * 100
        gen `topgrp'_sh_ccorw = ((ccorw_`topgrp' + forbes400_ccorw) / hweal_`topgrp') * 100
    }

    gen spec = "forbes3016replace"

    keep year spec *_sh_* 
    gen speclab = "11. Spec \#1 with Replace Top 400"
    append using `outdata'
    save `outdata', replace

    use `outdata', clear

end/*%>*/

capture program drop setup_suppC/*%<*/
program define setup_suppC

    tempfile outdata
    
    ****************************************************************************
    * Panel C: Growth in Top Wealth Shares through 2016 (p.p.)
    ****************************************************************************
    ****************************************************************************
    * 1. Baseline
    ****************************************************************************
    load_analysis_data szz

    #delimit ;
    keep if inlist(group, "All", "P99-100", "P99.9-100", "P99.99-100", "P99.999-100");

    replace group = cond(group == "P99.999-100", "_top0001",
                    cond(group == "P99.99-100", "_top001", 
                    cond(group == "P99.9-100", "_top01", 
                    cond(group == "P99-100", "_top1", 
                    cond(group == "All", "_total", "")))));
    drop if group == "";
    #delimit cr
    
    keep year group hweal_pref
    rename hweal_pref hweal
    reshape wide hweal, i(year) j(group, string)

    foreach topgrp in top1 top01 top001 top0001 {
        gen `topgrp'_share = (hweal_`topgrp' / hweal_total) * 100
    }
    keep if inlist(year, 1966, 1978, 1989, 2001, 2016)
    sort year 
    assert _N == 5

    foreach topgrp in top1 top01 top001 top0001 {
        local i = 1
        foreach year in 1966 1978 1989 2001 {
            gen d`topgrp'_`year' = `topgrp'_share[5] - `topgrp'_share[`i']
            local i = `i' + 1
        }
    }

    keep d*
    keep if _n == 1

    gen spec = "baseline"
    gen speclab = "1. Baseline"

    save `outdata'
    ****************************************************************************
    * 2. Equal returns
    ****************************************************************************
    load_analysis_data sz_v3

    #delimit ;
    keep if inlist(group, "All", "P99-100", "P99.9-100", "P99.99-100", "P99.999-100");

    replace group = cond(group == "P99.999-100", "_top0001",
                    cond(group == "P99.99-100", "_top001", 
                    cond(group == "P99.9-100", "_top01", 
                    cond(group == "P99-100", "_top1", 
                    cond(group == "All", "_total", "")))));
    drop if group == "";
    #delimit cr
    
    keep year group hweal_szv3
    rename hweal_szv3 hweal
    reshape wide hweal, i(year) j(group, string)

    foreach topgrp in top1 top01 top001 top0001 {
        gen `topgrp'_share = (hweal_`topgrp' / hweal_total) * 100
    }
    keep if inlist(year, 1966, 1978, 1989, 2001, 2016)
    sort year 
    assert _N == 5

    foreach topgrp in top1 top01 top001 top0001 {
        local i = 1
        foreach year in 1966 1978 1989 2001 {
            gen d`topgrp'_`year' = `topgrp'_share[5] - `topgrp'_share[`i']
            local i = `i' + 1
        }
    }

    keep d*
    keep if _n == 1

    gen spec = "equalrets"
    gen speclab = "2. Equal Returns"

    append using `outdata'
    save `outdata', replace
    ****************************************************************************
    * 3. SZ 2020
    ****************************************************************************
    load_analysis_data sz2020_top
    rename top*_pszes20 top*_share
    drop top10_share

    keep if inlist(year, 1966, 1978, 1989, 2001, 2016)
    sort year 
    assert _N == 5

    foreach topgrp in top1 top01 top001 top0001 {
        local i = 1
        foreach year in 1966 1978 1989 2001 {
            gen d`topgrp'_`year' = `topgrp'_share[5] - `topgrp'_share[`i']
            local i = `i' + 1
        }
    }

    keep d*
    keep if _n == 1

    gen spec = "revisedsz"
    gen speclab = "3. Revised Saez Zucman (2020)"

    append using `outdata'
    save `outdata', replace
    ****************************************************************************
    * 4. PSZ 2018
    ****************************************************************************
    load_analysis_data sz

    #delimit ;
    keep if inlist(group, "All", "P99-100", "P99.9-100", "P99.99-100", "P99.999-100");

    replace group = cond(group == "P99.999-100", "_top0001",
                    cond(group == "P99.99-100", "_top001", 
                    cond(group == "P99.9-100", "_top01", 
                    cond(group == "P99-100", "_top1", 
                    cond(group == "All", "_total", "")))));
    drop if group == "";
    #delimit cr
    
    keep year group hweal$sz_wlth_defn
    rename hweal$sz_wlth_defn hweal
    reshape wide hweal, i(year) j(group, string)

    foreach topgrp in top1 top01 top001 top0001 {
        gen `topgrp'_share = (hweal_`topgrp' / hweal_total) * 100
    }
    keep if inlist(year, 1966, 1978, 1989, 2001, 2016)
    sort year 
    assert _N == 5

    foreach topgrp in top1 top01 top001 top0001 {
        local i = 1
        foreach year in 1966 1978 1989 2001 {
            gen d`topgrp'_`year' = `topgrp'_share[5] - `topgrp'_share[`i']
            local i = `i' + 1
        }
    }

    keep d*
    keep if _n == 1

    gen spec = "psz"
    gen speclab = "4. Piketty Saez Zucman (2018)"

    append using `outdata'
    save `outdata', replace
    ****************************************************************************
    * 5, 6, 7, 8, 9: Alternative Aggregates
    ****************************************************************************
    * 5: 1 with unscaled pass-through
    * 6: 5 with missing pass-through
    * 7: 6 with unfunded pensions
    * 8: 7 with debt adjustments
    * 9: 8 with no student debt
    * 14: 2 with bond split (robustness)
    * 15: 1 with trimmed interest rates
    * 16: 8 with trimmed interest rates
    ****************************************************************************
    foreach rankspec in 3021 3022 3023 3024 3025 3018 3026 3027 {
        load_analysis_data altspec `rankspec' 

        #delimit ;
        keep if inlist(group, "All", "P99-100", "P99.9-100", "P99.99-100", "P99.999-100");

        replace group = cond(group == "P99.999-100", "_top0001",
                        cond(group == "P99.99-100", "_top001", 
                        cond(group == "P99.9-100", "_top01", 
                        cond(group == "P99-100", "_top1", 
                        cond(group == "All", "_total", "")))));
        drop if group == "";
        #delimit cr
        
        keep year group hweal_alt
        rename hweal_alt hweal
        reshape wide hweal, i(year) j(group, string)

        foreach topgrp in top1 top01 top001 top0001 {
            gen `topgrp'_share = (hweal_`topgrp' / hweal_total) * 100
        }
        keep if inlist(year, 1966, 1978, 1989, 2001, 2016)
        sort year 
        assert _N == 5

        foreach topgrp in top1 top01 top001 top0001 {
            local i = 1
            foreach year in 1966 1978 1989 2001 {
                gen d`topgrp'_`year' = `topgrp'_share[5] - `topgrp'_share[`i']
                local i = `i' + 1
            }
        }

        keep d*
        keep if _n == 1

        gen spec = "alt`rankspec'"

        if (`rankspec' == 3021) {
            gen speclab = "5. Spec \#1 with Unscaled Pass-through"
        }
        else if (`rankspec' == 3022) {
            gen speclab = "6. Spec \#5 with Missing Pass-through"
        }
        else if (`rankspec' == 3023) {
            gen speclab = "7. Spec \#6 with Unfunded Pensions"
        }
        else if (`rankspec' == 3024) {
            gen speclab = "8. Spec \#7 with Debt Adjustments"
        }
        else if (`rankspec' == 3025) {
            gen speclab = "9. Spec \#8 with No Student Debt"
        }
        else if (`rankspec' == 3018) {
            gen speclab = "14. Spec \#2 with Bond Split"
        }
        else if (`rankspec' == 3026) {
            gen speclab = "15. Spec \#1 with Low Boutique"
        }
        else if (`rankspec' == 3027) {
            gen speclab = "16. Spec \#8 with Low Boutique"
        }
        append using `outdata'
        save `outdata', replace
    }

    ****************************************************************************
    * 10, 11, 12, 13, 17: Forbes specs
    ****************************************************************************
    * 10: Baseline with BHV
    * 11: Baseline with Replace
    * 12: 6 with BHV
    * 13: 8 with BHV
    * 17: 13 with BHV
    ****************************************************************************
    foreach rankspec in 3016 3022 3024 3027 {
        load_analysis_data szz_forbes `rankspec'

        foreach topgrp in top1 top01 top001 top0001 {
            gen `topgrp'_share = `topgrp'share_adjusted 
        }

        keep if inlist(year, 1982, 1989, 2001, 2016)
        keep year *_share 
        sort year 
        assert _N == 4

        foreach topgrp in top1 top01 top001 top0001 {
            local i = 1
            foreach year in 1982 1989 2001 {
                gen d`topgrp'_`year' = `topgrp'_share[4] - `topgrp'_share[`i']
                local i = `i' + 1
            }
        }

        keep d*
        keep if _n == 1

        gen spec = "forbes`rankspec'"

        if ("`rankspec'" == "3016") {
            gen speclab = "10. Spec \#1 with Reweighting"
        }
        else if ("`rankspec'" == "3022") {
            gen speclab = "12. Spec \#6 with Reweighting"
        }
        else if ("`rankspec'" == "3024") {
            gen speclab = "13. Spec \#8 with Reweighting"
        }
        else if ("`rankspec'" == "3027") {
            gen speclab = "17. Spec \#13 with Low Boutique"
        }
        append using `outdata'
        save `outdata', replace
    }

    load_analysis_data szz_forbes 3016

    foreach topgrp in top1 top01 top001 top0001 {
        gen `topgrp'_share = `topgrp'share_replaced 
    }

    keep if inlist(year, 1982, 1989, 2001, 2016)
    keep year *_share 
    sort year 
    assert _N == 4

    foreach topgrp in top1 top01 top001 top0001 {
        local i = 1
        foreach year in 1982 1989 2001 {
            gen d`topgrp'_`year' = `topgrp'_share[4] - `topgrp'_share[`i']
            local i = `i' + 1
        }
    }

    keep d*
    keep if _n == 1

    gen spec = "forbes3016replace"
    gen speclab = "11. Spec \#1 with Replace Top 400"
    append using `outdata'
    save `outdata', replace

    use `outdata', clear

end/*%>*/

capture program drop setup_suppD/*%<*/
program define setup_suppD

    tempfile outdata
    
    ****************************************************************************
    * Panel D: Growth in Top Portfolio Shares in 2001--2016 (%)
    ****************************************************************************
    ****************************************************************************
    * 1. Baseline
    ****************************************************************************
    load_taxdata, rankspec($preferred_defn_late) startyr(2001) endyr(2016)

    #delimit ;
    keep if inlist(group, "All", "P99-100", "P99.9-100", "P99.99-100", "P99.999-100");

    replace group = cond(group == "P99.999-100", "_top0001",
                    cond(group == "P99.99-100", "_top001", 
                    cond(group == "P99.9-100", "_top01", 
                    cond(group == "P99-100", "_top1", 
                    cond(group == "All", "_total", "")))));
    drop if group == "";
    #delimit cr
    
    keep year group hweal$preferred_defn_late hwfix_base pthru_base ccorw_base 
    rename hweal$preferred_defn_late hweal
    rename *_base *
    reshape wide hweal hwfix pthru ccorw, i(year) j(group, string)

    foreach topgrp in top1 top01 top001 top0001 {
        gen `topgrp'_sh_hwfix = (hwfix_`topgrp' / hweal_`topgrp') * 100
        gen `topgrp'_sh_pthru = (pthru_`topgrp' / hweal_`topgrp') * 100
        gen `topgrp'_sh_ccorw = (ccorw_`topgrp' / hweal_`topgrp') * 100
    }
    keep if inlist(year, 2001, 2016)
    sort year
    assert _N == 2

    foreach topgrp in top1 top01 top001 top0001 {
        local i = 1
        gen d`topgrp'_sh_hwfix = `topgrp'_sh_hwfix[2] - `topgrp'_sh_hwfix[1]
        gen d`topgrp'_sh_pthru = `topgrp'_sh_pthru[2] - `topgrp'_sh_pthru[1]
        gen d`topgrp'_sh_ccorw = `topgrp'_sh_ccorw[2] - `topgrp'_sh_ccorw[1]
        local i = `i' + 1
    }

    keep d*
    keep if _n == 1

    gen spec = "baseline"
    gen speclab = "1. Baseline"

    save `outdata'
    ****************************************************************************
    * 2. Equal returns
    ****************************************************************************
    load_taxdata, rankspec($sz_wlth_defn_late_v3) startyr(2001) endyr(2016)

    #delimit ;
    keep if inlist(group, "All", "P99-100", "P99.9-100", "P99.99-100", "P99.999-100");

    replace group = cond(group == "P99.999-100", "_top0001",
                    cond(group == "P99.99-100", "_top001", 
                    cond(group == "P99.9-100", "_top01", 
                    cond(group == "P99-100", "_top1", 
                    cond(group == "All", "_total", "")))));
    drop if group == "";
    #delimit cr
    
    keep year group hweal$sz_wlth_defn_late_v3 hwfix_equ pthru_equ ccorw_equ
    rename hweal$sz_wlth_defn_late_v3 hweal
    rename *_equ *
    reshape wide hweal hwfix pthru ccorw, i(year) j(group, string)

    foreach topgrp in top1 top01 top001 top0001 {
        gen `topgrp'_sh_hwfix = (hwfix_`topgrp' / hweal_`topgrp') * 100
        gen `topgrp'_sh_pthru = (pthru_`topgrp' / hweal_`topgrp') * 100
        gen `topgrp'_sh_ccorw = (ccorw_`topgrp' / hweal_`topgrp') * 100
    }
    keep if inlist(year, 2001, 2016)
    sort year
    assert _N == 2

    foreach topgrp in top1 top01 top001 top0001 {
        local i = 1
        gen d`topgrp'_sh_hwfix = `topgrp'_sh_hwfix[2] - `topgrp'_sh_hwfix[1]
        gen d`topgrp'_sh_pthru = `topgrp'_sh_pthru[2] - `topgrp'_sh_pthru[1]
        gen d`topgrp'_sh_ccorw = `topgrp'_sh_ccorw[2] - `topgrp'_sh_ccorw[1]
        local i = `i' + 1
    }

    keep d*
    keep if _n == 1

    gen spec = "equalrets"
    gen speclab = "2. Equal Returns"

    append using `outdata'
    save `outdata', replace
    ****************************************************************************
    * 3. SZ 2020
    ****************************************************************************
    load_analysis_data sz2020_port
    keep if year == 2001 | year == 2016

    foreach topgrp in top1 top01 top001 {
        gen `topgrp'_sh_hwfix = (fix_`topgrp' / wlth_`topgrp') * 100
        gen `topgrp'_sh_pthru = .
        gen `topgrp'_sh_ccorw = .
    }
    sort year
    assert _N == 2

    foreach topgrp in top1 top01 top001 {
        local i = 1
        gen d`topgrp'_sh_hwfix = `topgrp'_sh_hwfix[2] - `topgrp'_sh_hwfix[1]
        gen d`topgrp'_sh_pthru = `topgrp'_sh_pthru[2] - `topgrp'_sh_pthru[1]
        gen d`topgrp'_sh_ccorw = `topgrp'_sh_ccorw[2] - `topgrp'_sh_ccorw[1]
        local i = `i' + 1
    }

    keep d*
    keep if _n == 1

    gen spec = "revisedsz"
    gen speclab = "3. Revised Saez Zucman (2020)"

    append using `outdata'
    save `outdata', replace
    ****************************************************************************
    * 4. PSZ 2018
    ****************************************************************************
    load_taxdata, rankspec($sz_wlth_defn) startyr(2001) endyr(2016)

    #delimit ;
    keep if inlist(group, "All", "P99-100", "P99.9-100", "P99.99-100", "P99.999-100");

    replace group = cond(group == "P99.999-100", "_top0001",
                    cond(group == "P99.99-100", "_top001", 
                    cond(group == "P99.9-100", "_top01", 
                    cond(group == "P99-100", "_top1", 
                    cond(group == "All", "_total", "")))));
    drop if group == "";
    #delimit cr
    
    keep year group hweal$sz_wlth_defn hwfix_equ pthru_equ ccorw_equ
    rename hweal$sz_wlth_defn hweal
    rename *_equ *
    reshape wide hweal hwfix pthru ccorw, i(year) j(group, string)

    foreach topgrp in top1 top01 top001 top0001 {
        gen `topgrp'_sh_hwfix = (hwfix_`topgrp' / hweal_`topgrp') * 100
        gen `topgrp'_sh_pthru = (pthru_`topgrp' / hweal_`topgrp') * 100
        gen `topgrp'_sh_ccorw = (ccorw_`topgrp' / hweal_`topgrp') * 100
    }
    keep if inlist(year, 2001, 2016)
    sort year
    assert _N == 2

    foreach topgrp in top1 top01 top001 top0001 {
        local i = 1
        gen d`topgrp'_sh_hwfix = `topgrp'_sh_hwfix[2] - `topgrp'_sh_hwfix[1]
        gen d`topgrp'_sh_pthru = `topgrp'_sh_pthru[2] - `topgrp'_sh_pthru[1]
        gen d`topgrp'_sh_ccorw = `topgrp'_sh_ccorw[2] - `topgrp'_sh_ccorw[1]
        local i = `i' + 1
    }

    keep d*
    keep if _n == 1

    gen spec = "psz"
    gen speclab = "4. Piketty Saez Zucman (2018)"

    append using `outdata'
    save `outdata', replace
    ****************************************************************************
    * 5, 6, 7, 8, 9: Alternative Aggregates
    ****************************************************************************
    * 5: 1 with unscaled pass-through
    * 6: 5 with missing pass-through
    * 7: 6 with unfunded pensions
    * 8: 7 with debt adjustments
    * 9: 8 with no student debt
    * 14: 2 with bond split (robustness)
    * 15: 1 with trimmed interest rates
    * 16: 8 with trimmed interest rates
    ****************************************************************************
    foreach rankspec in 3021 3022 3023 3024 3025 3018 3026 3027 {
        load_taxdata, rankspec(`rankspec') startyr(2001) endyr(2016)

        #delimit ;
        keep if inlist(group, "All", "P99-100", "P99.9-100", "P99.99-100", "P99.999-100");

        replace group = cond(group == "P99.999-100", "_top0001",
                        cond(group == "P99.99-100", "_top001", 
                        cond(group == "P99.9-100", "_top01", 
                        cond(group == "P99-100", "_top1", 
                        cond(group == "All", "_total", "")))));
        drop if group == "";
        #delimit cr
        
        * Note hwfix_base = hwfix_pref and ccorw_base = ccorw_pref
        if (`rankspec' == 3021) {
            keep year group hweal`rankspec' hwfix_pref pthru_pref ccorw_pref ///
                missing_scorp missing_pship
            replace pthru_pref = pthru_pref - missing_scorp - missing_pship
            rename *_pref *
            drop missing_scorp missing_pship
        }
        else if inlist(`rankspec', 3022, 3023, 3024, 3025) {
            keep year group hweal`rankspec' hwfix_pref pthru_pref ccorw_pref 
            rename *_pref *
        }
		else if inlist(`rankspec', 3026) {
            keep year group hweal`rankspec' hwfix_pref pthru_base ccorw_pref ///
                taxbond_info_trim taxbond_info
            replace hwfix_pref = hwfix_pref - taxbond_info + taxbond_info_trim
            rename *_pref *
            drop taxbond*
			rename *_base *
        }
        else if inlist(`rankspec', 3027) {
            keep year group hweal`rankspec' hwfix_pref pthru_pref ccorw_pref ///
                taxbond_info_trim taxbond_info
            replace hwfix_pref = hwfix_pref - taxbond_info + taxbond_info_trim
            rename *_pref *
            drop taxbond*
        }
        else if inlist(`rankspec', 3018) {
            keep year group hweal`rankspec' hwfix_equal pthru_equ ccorw_equal 
            rename *_equal *
            rename *_equ *
        }
        rename hweal`rankspec' hweal
        reshape wide hweal hwfix pthru ccorw, i(year) j(group, string)

        foreach topgrp in top1 top01 top001 top0001 {
            gen `topgrp'_sh_hwfix = (hwfix_`topgrp' / hweal_`topgrp') * 100
            gen `topgrp'_sh_pthru = (pthru_`topgrp' / hweal_`topgrp') * 100
            gen `topgrp'_sh_ccorw = (ccorw_`topgrp' / hweal_`topgrp') * 100
        }
        keep if inlist(year, 2001, 2016)
        sort year
        assert _N == 2

        foreach topgrp in top1 top01 top001 top0001 {
            local i = 1
            gen d`topgrp'_sh_hwfix = `topgrp'_sh_hwfix[2] - `topgrp'_sh_hwfix[1]
            gen d`topgrp'_sh_pthru = `topgrp'_sh_pthru[2] - `topgrp'_sh_pthru[1]
            gen d`topgrp'_sh_ccorw = `topgrp'_sh_ccorw[2] - `topgrp'_sh_ccorw[1]
            local i = `i' + 1
        }

        keep d*
        keep if _n == 1

        gen spec = "alt`rankspec'"

        if (`rankspec' == 3021) {
            gen speclab = "5. Spec \#1 with Unscaled Pass-through"
        }
        else if (`rankspec' == 3022) {
            gen speclab = "6. Spec \#5 with Missing Pass-through"
        }
        else if (`rankspec' == 3023) {
            gen speclab = "7. Spec \#6 with Unfunded Pensions"
        }
        else if (`rankspec' == 3024) {
            gen speclab = "8. Spec \#7 with Debt Adjustments"
        }
        else if (`rankspec' == 3025) {
            gen speclab = "9. Spec \#8 with No Student Debt"
        }
        else if (`rankspec' == 3018) {
            gen speclab = "14. Spec \#2 with Bond Split"
        }
        else if (`rankspec' == 3026) {
            gen speclab = "15. Spec \#1 with Low Boutique"
        }
        else if (`rankspec' == 3027) {
            gen speclab = "16. Spec \#8 with Low Boutique"
        }
        append using `outdata'
        save `outdata', replace
    }
    ****************************************************************************
    * 10, 11, 12, 13: Forbes specs
    ****************************************************************************
    * 10: Baseline with BHV
    * 11: Baseline with Replace
    * 12: 6 with BHV
    * 13: 8 with BHV
    * 17: 13 with BHV
    ****************************************************************************
    foreach rankspec in 3016 3022 3024 3027 {
        load_analysis_data szz_forbes `rankspec'

        if inlist(`rankspec', 3022, 3024, 3027) {
            keep year hweal_* hwfix_pref_*_adjusted pthru_pref_*_adjusted ccorw_pref_*_adjusted ///
                forbes400_hwfix_adjusted forbes400_pthru_adjusted forbes400_ccorw_adjusted
            rename *_adjusted *
            rename hwfix_pref_* hwfix_*
            rename pthru_pref_* pthru_*
            rename ccorw_pref_* ccorw_*
        }
        else if inlist(`rankspec', 3016) {
            keep year hweal_* hwfix_base_*_adjusted pthru_base_*_adjusted ccorw_base_*_adjusted ///
                forbes400_hwfix_adjusted forbes400_pthru_adjusted forbes400_ccorw_adjusted
            rename *_adjusted *
            rename hwfix_base_* hwfix_*
            rename pthru_base_* pthru_*
            rename ccorw_base_* ccorw_*
        }

        foreach topgrp in top1 top01 top001 top0001 {
            gen `topgrp'_sh_hwfix = ((hwfix_`topgrp' + forbes400_hwfix) / hweal_`topgrp') * 100
            gen `topgrp'_sh_pthru = ((pthru_`topgrp' + forbes400_pthru) / hweal_`topgrp') * 100
            gen `topgrp'_sh_ccorw = ((ccorw_`topgrp' + forbes400_ccorw) / hweal_`topgrp') * 100
        }
        keep if inlist(year, 2001, 2016)
        sort year 
        assert _N == 2

        foreach topgrp in top1 top01 top001 top0001 {
            local i = 1
            gen d`topgrp'_sh_hwfix = `topgrp'_sh_hwfix[2] - `topgrp'_sh_hwfix[1]
            gen d`topgrp'_sh_pthru = `topgrp'_sh_pthru[2] - `topgrp'_sh_pthru[1]
            gen d`topgrp'_sh_ccorw = `topgrp'_sh_ccorw[2] - `topgrp'_sh_ccorw[1]
            local i = `i' + 1
        }

        keep d*
        keep if _n == 1

        gen spec = "forbes`rankspec'"

        keep spec *_sh_* 
        if ("`rankspec'" == "3016") {
            gen speclab = "10. Spec \#1 with Reweighting"
        }
        else if ("`rankspec'" == "3022") {
            gen speclab = "12. Spec \#6 with Reweighting"
        }
        else if ("`rankspec'" == "3024") {
            gen speclab = "13. Spec \#8 with Reweighting"
        }
        else if ("`rankspec'" == "3027") {
            gen speclab = "17. Spec \#13 with Low Boutique"
        }
        append using `outdata'
        save `outdata', replace
    }

    load_analysis_data szz_forbes 3016

	keep year hweal_* hwfix_base_*_replaced pthru_base_*_replaced ccorw_base_*_replaced ///
        forbes400_hwfix_replaced forbes400_pthru_replaced forbes400_ccorw_replaced
    rename *_replaced *
    rename hwfix_base_* hwfix_*
    rename pthru_base_* pthru_*
    rename ccorw_base_* ccorw_*

    foreach topgrp in top1 top01 top001 top0001 {
        gen `topgrp'_sh_hwfix = ((hwfix_`topgrp' + forbes400_hwfix) / hweal_`topgrp') * 100
        gen `topgrp'_sh_pthru = ((pthru_`topgrp' + forbes400_pthru) / hweal_`topgrp') * 100
        gen `topgrp'_sh_ccorw = ((ccorw_`topgrp' + forbes400_ccorw) / hweal_`topgrp') * 100
    }
    keep if inlist(year, 2001, 2016)
    sort year 
    assert _N == 2

    foreach topgrp in top1 top01 top001 top0001 {
        local i = 1
        gen d`topgrp'_sh_hwfix = `topgrp'_sh_hwfix[2] - `topgrp'_sh_hwfix[1]
        gen d`topgrp'_sh_pthru = `topgrp'_sh_pthru[2] - `topgrp'_sh_pthru[1]
        gen d`topgrp'_sh_ccorw = `topgrp'_sh_ccorw[2] - `topgrp'_sh_ccorw[1]
        local i = `i' + 1
    }

    keep d*
    keep if _n == 1

    gen spec = "forbes3016replace"

    keep spec *_sh_* 
    gen speclab = "11. Spec \#1 with Replace Top 400"
    append using `outdata'
    save `outdata', replace

    use `outdata', clear

end/*%>*/

capture program drop table_supplemental
program define table_supplemental

    tempfile suppA
    setup_suppA
    save `suppA'

    tempfile suppB
    setup_suppB
    save `suppB'

    tempfile suppC
    setup_suppC
    save `suppC'

    tempfile suppD
    setup_suppD
    save `suppD'

    ****************************************************************************
    * Make Table Fragments (Panel A)%<
    ****************************************************************************
    clear
    set obs 1

    gen tab = "\begin{tabular}{lcccccccccc}"
    gen toprule = "\toprule"
    gen midrule = "\midrule"
    gen botrule = "\bottomrule"
    gen cmidrule = "\cmidrule(lr){2-3} \cmidrule(lr){4-7} \cmidrule(lr){8-11}"
    gen vsp = "\vspace{.03in}" in 1
    gen end = "\end{tabular}" in 1

    gen c = " "

    gen title1 = "\multicolumn{11}{l}{Panel A. Top Wealth Shares and Per-Capita Wealth in 2016}"
    gen title2 = "& \multicolumn{2}{c}{Population} & \multicolumn{4}{c}{Top Wealth Share (\%)} &   \multicolumn{4}{c}{Per-Capita Wealth ({\\$}M)}"
    gen title3 = "& Agg ({\\$}T) & Per Cap ({\\$}K) & 1\% & 0.1\% & 0.01\% & 0.001\% & 1\% & 0.1\% & 0.01\% & 0.001\%"

    gen sec1 = "\textbf{Benchmark Series} & & & & & & & & & &"
    gen sec2 = "\textbf{Alternative Aggregates} & & & & & & & & & &"
    gen sec3 = "\textbf{Forbes Augmentation} & & & & & & & & & &"
    gen sec4 = "\textbf{Robustness Series} & & & & & & & & & &"

    append using `suppA'
    gen total_percap = hweal_total / n_total

    local tabname = "supp_series_panelA"
    
    * Top matter
    listtex tab if _n == 1 using "`tabname'.tex", replace rstyle(none)
    listtex title1 if _n == 1, appendto("`tabname'.tex") rstyle(tabular)
    listtex toprule if _n == 1, appendto("`tabname'.tex") rstyle(none)
    listtex title2 if _n == 1, appendto("`tabname'.tex") rstyle(tabular)
    listtex cmidrule if _n == 1, appendto("`tabname'.tex") rstyle(none)
    listtex title3 if _n == 1, appendto("`tabname'.tex") rstyle(tabular)
    listtex midrule if _n==1 , appendto("`tabname'.tex") rstyle(none)    

    * Formatting
    replace hweal_total = round(hweal_total * 1e-12, 0.1)
    foreach v of varlist top*_share {
        replace `v' = round(`v', 0.1)
        format `v' %9.1f
    }
    foreach v of varlist top*_percap {
        replace `v' = round(`v' * 1e-6, 0.1)
        format `v' %9.1fc
    }
    replace total_percap = round(total_percap * 1e-3, 0.1)
    format total_percap %9.1f

    * Sections
    #delimit ;
    listtex sec1 if _n == 1, appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab hweal_total total_percap top1_share top01_share top001_share top0001_share
        top1_percap top01_percap top001_percap top0001_percap 
        if spec == "baseline", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab hweal_total total_percap top1_share top01_share top001_share top0001_share
        top1_percap top01_percap top001_percap top0001_percap 
        if spec == "equalrets", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab hweal_total total_percap top1_share top01_share top001_share top0001_share
        top1_percap top01_percap top001_percap top0001_percap 
        if spec == "revisedsz", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab hweal_total total_percap top1_share top01_share top001_share top0001_share
        top1_percap top01_percap top001_percap top0001_percap 
        if spec == "psz", appendto("`tabname'.tex") rstyle(tabular);
    listtex vsp if _n == 1, appendto("`tabname'.tex") rstyle(tabular);

    listtex sec2 if _n == 1, appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab hweal_total total_percap top1_share top01_share top001_share top0001_share
        top1_percap top01_percap top001_percap top0001_percap 
        if spec == "alt3021", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab hweal_total total_percap top1_share top01_share top001_share top0001_share
        top1_percap top01_percap top001_percap top0001_percap 
        if spec == "alt3022", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab hweal_total total_percap top1_share top01_share top001_share top0001_share
        top1_percap top01_percap top001_percap top0001_percap 
        if spec == "alt3023", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab hweal_total total_percap top1_share top01_share top001_share top0001_share
        top1_percap top01_percap top001_percap top0001_percap 
        if spec == "alt3024", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab hweal_total total_percap top1_share top01_share top001_share top0001_share
        top1_percap top01_percap top001_percap top0001_percap 
        if spec == "alt3025", appendto("`tabname'.tex") rstyle(tabular);
    listtex vsp if _n == 1, appendto("`tabname'.tex") rstyle(tabular);

    listtex sec3 if _n == 1, appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab hweal_total total_percap top1_share top01_share top001_share top0001_share
        top1_percap top01_percap top001_percap top0001_percap 
        if spec == "forbes3016", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab hweal_total total_percap top1_share top01_share top001_share top0001_share
        top1_percap top01_percap top001_percap top0001_percap 
        if spec == "forbes3016replace", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab hweal_total total_percap top1_share top01_share top001_share top0001_share
        top1_percap top01_percap top001_percap top0001_percap 
        if spec == "forbes3022", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab hweal_total total_percap top1_share top01_share top001_share top0001_share
        top1_percap top01_percap top001_percap top0001_percap 
        if spec == "forbes3024", appendto("`tabname'.tex") rstyle(tabular);
    listtex vsp if _n == 1, appendto("`tabname'.tex") rstyle(tabular);

    listtex sec4 if _n == 1, appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab hweal_total total_percap top1_share top01_share top001_share top0001_share
        top1_percap top01_percap top001_percap top0001_percap 
        if spec == "alt3018", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab hweal_total total_percap top1_share top01_share top001_share top0001_share
        top1_percap top01_percap top001_percap top0001_percap 
        if spec == "alt3026", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab hweal_total total_percap top1_share top01_share top001_share top0001_share
        top1_percap top01_percap top001_percap top0001_percap 
        if spec == "alt3027", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab hweal_total total_percap top1_share top01_share top001_share top0001_share
        top1_percap top01_percap top001_percap top0001_percap 
        if spec == "forbes3027", appendto("`tabname'.tex") rstyle(tabular);
    #delimit cr

    qui listtex botrule if _n == 1, appendto("`tabname'.tex") rstyle(none)  
    qui listtex end if _n == 1, appendto("`tabname'.tex") rstyle(none)/*%>*/
    ****************************************************************************

    ****************************************************************************
    * Make Table Fragments (Panel B)%<
    ****************************************************************************
    clear
    set obs 1

    gen tab = "\begin{tabular}{lccccccccc}"
    gen toprule = "\toprule"
    gen midrule = "\midrule"
    gen botrule = "\bottomrule"
    gen cmidrule = "\cmidrule(lr){2-4} \cmidrule(lr){5-7} \cmidrule(lr){8-10}"
    gen vsp = "\vspace{.03in}" in 1
    gen end = "\end{tabular}" in 1

    gen c = " "

    gen title1 = "\multicolumn{10}{l}{Panel B. Top Portfolio Shares in 2016 (\%)}"
    gen title2 = "& \multicolumn{3}{c}{Top 0.1\%} &   \multicolumn{3}{c}{Top 0.01\%} &   \multicolumn{3}{c}{Top 0.001\%}"
    gen title3 = "& Fix & Pthru & C-corp & Fix & Pthru & C-corp & Fix & Pthru & C-corp"

    gen sec1 = "\textbf{Benchmark Series} & & & & & & & & &"
    gen sec2 = "\textbf{Alternative Aggregates} & & & & & & & & &"
    gen sec3 = "\textbf{Forbes Augmentation} & & & & & & & & &"
    gen sec4 = "\textbf{Robustness Series} & & & & & & & & &"

    append using `suppB'

    local tabname = "supp_series_panelB"
    
    * Top matter
    listtex tab if _n == 1 using "`tabname'.tex", replace rstyle(none)
    listtex title1 if _n == 1, appendto("`tabname'.tex") rstyle(tabular)
    listtex toprule if _n == 1, appendto("`tabname'.tex") rstyle(none)
    listtex title2 if _n == 1, appendto("`tabname'.tex") rstyle(tabular)
    listtex cmidrule if _n == 1, appendto("`tabname'.tex") rstyle(none)
    listtex title3 if _n == 1, appendto("`tabname'.tex") rstyle(tabular)
    listtex midrule if _n==1 , appendto("`tabname'.tex") rstyle(none)    

    * Formatting
    foreach v of varlist top*_sh* {
        replace `v' = round(`v', 0.1)
        format `v' %3.1f
    }

    * Sections
    #delimit ;
    listtex sec1 if _n == 1, appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab top01_sh_hwfix top01_sh_pthru top01_sh_ccorw top001_sh_hwfix top001_sh_pthru top001_sh_ccorw
        top0001_sh_hwfix top0001_sh_pthru top0001_sh_ccorw 
        if spec == "baseline", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab top01_sh_hwfix top01_sh_pthru top01_sh_ccorw top001_sh_hwfix top001_sh_pthru top001_sh_ccorw
        top0001_sh_hwfix top0001_sh_pthru top0001_sh_ccorw 
        if spec == "equalrets", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab top01_sh_hwfix top01_sh_pthru top01_sh_ccorw top001_sh_hwfix top001_sh_pthru top001_sh_ccorw
        top0001_sh_hwfix top0001_sh_pthru top0001_sh_ccorw 
        if spec == "revisedsz", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab top01_sh_hwfix top01_sh_pthru top01_sh_ccorw top001_sh_hwfix top001_sh_pthru top001_sh_ccorw
        top0001_sh_hwfix top0001_sh_pthru top0001_sh_ccorw 
        if spec == "psz", appendto("`tabname'.tex") rstyle(tabular);
    listtex vsp if _n == 1, appendto("`tabname'.tex") rstyle(tabular);

    listtex sec2 if _n == 1, appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab top01_sh_hwfix top01_sh_pthru top01_sh_ccorw top001_sh_hwfix top001_sh_pthru top001_sh_ccorw
        top0001_sh_hwfix top0001_sh_pthru top0001_sh_ccorw 
        if spec == "alt3021", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab top01_sh_hwfix top01_sh_pthru top01_sh_ccorw top001_sh_hwfix top001_sh_pthru top001_sh_ccorw
        top0001_sh_hwfix top0001_sh_pthru top0001_sh_ccorw 
        if spec == "alt3022", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab top01_sh_hwfix top01_sh_pthru top01_sh_ccorw top001_sh_hwfix top001_sh_pthru top001_sh_ccorw
        top0001_sh_hwfix top0001_sh_pthru top0001_sh_ccorw 
        if spec == "alt3023", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab top01_sh_hwfix top01_sh_pthru top01_sh_ccorw top001_sh_hwfix top001_sh_pthru top001_sh_ccorw
        top0001_sh_hwfix top0001_sh_pthru top0001_sh_ccorw 
        if spec == "alt3024", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab top01_sh_hwfix top01_sh_pthru top01_sh_ccorw top001_sh_hwfix top001_sh_pthru top001_sh_ccorw
        top0001_sh_hwfix top0001_sh_pthru top0001_sh_ccorw 
        if spec == "alt3025", appendto("`tabname'.tex") rstyle(tabular);
    listtex vsp if _n == 1, appendto("`tabname'.tex") rstyle(tabular);

    listtex sec3 if _n == 1, appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab 
        top01_sh_hwfix top01_sh_pthru top01_sh_ccorw top001_sh_hwfix top001_sh_pthru top001_sh_ccorw
        top0001_sh_hwfix top0001_sh_pthru top0001_sh_ccorw
        if spec == "forbes3016", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab 
        top01_sh_hwfix top01_sh_pthru top01_sh_ccorw top001_sh_hwfix top001_sh_pthru top001_sh_ccorw
        top0001_sh_hwfix top0001_sh_pthru top0001_sh_ccorw
        if spec == "forbes3016replace", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab 
        top01_sh_hwfix top01_sh_pthru top01_sh_ccorw top001_sh_hwfix top001_sh_pthru top001_sh_ccorw
        top0001_sh_hwfix top0001_sh_pthru top0001_sh_ccorw
        if spec == "forbes3022", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab 
        top01_sh_hwfix top01_sh_pthru top01_sh_ccorw top001_sh_hwfix top001_sh_pthru top001_sh_ccorw
        top0001_sh_hwfix top0001_sh_pthru top0001_sh_ccorw
        if spec == "forbes3024", appendto("`tabname'.tex") rstyle(tabular);
    listtex vsp if _n == 1, appendto("`tabname'.tex") rstyle(tabular);   

    listtex sec4 if _n == 1, appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab top01_sh_hwfix top01_sh_pthru top01_sh_ccorw top001_sh_hwfix top001_sh_pthru top001_sh_ccorw
        top0001_sh_hwfix top0001_sh_pthru top0001_sh_ccorw 
        if spec == "alt3018", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab top01_sh_hwfix top01_sh_pthru top01_sh_ccorw top001_sh_hwfix top001_sh_pthru top001_sh_ccorw
        top0001_sh_hwfix top0001_sh_pthru top0001_sh_ccorw 
        if spec == "alt3026", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab top01_sh_hwfix top01_sh_pthru top01_sh_ccorw top001_sh_hwfix top001_sh_pthru top001_sh_ccorw
        top0001_sh_hwfix top0001_sh_pthru top0001_sh_ccorw 
        if spec == "alt3027", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab top01_sh_hwfix top01_sh_pthru top01_sh_ccorw top001_sh_hwfix top001_sh_pthru top001_sh_ccorw
        top0001_sh_hwfix top0001_sh_pthru top0001_sh_ccorw 
        if spec == "forbes3027", appendto("`tabname'.tex") rstyle(tabular);
    #delimit cr

    qui listtex botrule if _n == 1, appendto("`tabname'.tex") rstyle(none)  
    qui listtex end if _n == 1, appendto("`tabname'.tex") rstyle(none)/*%>*/
    ****************************************************************************

    ****************************************************************************
    * Make Table Fragments (Panel C)%<
    ****************************************************************************
    clear
    set obs 1

    gen tab = "\begin{tabular}{lcccccccccccc}"
    gen toprule = "\toprule"
    gen midrule = "\midrule"
    gen botrule = "\bottomrule"
    gen cmidrule = "\cmidrule(lr){2-5} \cmidrule(lr){6-9} \cmidrule(lr){10-13}"
    gen vsp = "\vspace{.03in}" in 1
    gen end = "\end{tabular}" in 1

    gen c = " "

    gen title1 = "\multicolumn{11}{l}{Panel C. Growth in Top Wealth Shares through 2016 (p.p.)}"
    gen title2 = "& \multicolumn{4}{c}{Top 0.1\% Since} &   \multicolumn{4}{c}{Top 0.01\% Since} &   \multicolumn{4}{c}{Top 0.001\% Since}"
    gen title3 = "& 1966 & 1978 & 1989 & 2001 & 1966 & 1978 & 1989 & 2001 & 1966 & 1978 & 1989 & 2001"

    gen title3alt = "& 1966 & 1982 & 1989 & 2001 & 1966 & 1982 & 1989 & 2001 & 1966 & 1982 & 1989 & 2001"

    gen sec1 = "\textbf{Benchmark Series} & & & & & & & & & &"
    gen sec2 = "\textbf{Alternative Aggregates} & & & & & & & & & &"
    gen sec3 = "\textbf{Forbes Augmentation} & & & & & & & & & &"
    gen sec4 = "\textbf{Robustness Series} & & & & & & & & & &"

    append using `suppC'

    local tabname = "supp_series_panelC"
    
    * Top matter
    listtex tab if _n == 1 using "`tabname'.tex", replace rstyle(none)
    listtex title1 if _n == 1, appendto("`tabname'.tex") rstyle(tabular)
    listtex toprule if _n == 1, appendto("`tabname'.tex") rstyle(none)
    listtex title2 if _n == 1, appendto("`tabname'.tex") rstyle(tabular)
    listtex cmidrule if _n == 1, appendto("`tabname'.tex") rstyle(none)
    listtex title3 if _n == 1, appendto("`tabname'.tex") rstyle(tabular)
    listtex midrule if _n==1 , appendto("`tabname'.tex") rstyle(none)    

    * Formatting
    foreach v of varlist dtop* {
        replace `v' = round(`v', 0.1)
        format `v' %3.1f
    }

    * Sections
    #delimit ;
    listtex sec1 if _n == 1, appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab  
        dtop01_1966 dtop01_1978 dtop01_1989 dtop01_2001
        dtop001_1966 dtop001_1978 dtop001_1989 dtop001_2001
        dtop0001_1966 dtop0001_1978 dtop0001_1989 dtop0001_2001
        if spec == "baseline", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab  
        dtop01_1966 dtop01_1978 dtop01_1989 dtop01_2001
        dtop001_1966 dtop001_1978 dtop001_1989 dtop001_2001
        dtop0001_1966 dtop0001_1978 dtop0001_1989 dtop0001_2001
        if spec == "equalrets", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab  
        dtop01_1966 dtop01_1978 dtop01_1989 dtop01_2001
        dtop001_1966 dtop001_1978 dtop001_1989 dtop001_2001
        dtop0001_1966 dtop0001_1978 dtop0001_1989 dtop0001_2001
        if spec == "revisedsz", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab  
        dtop01_1966 dtop01_1978 dtop01_1989 dtop01_2001
        dtop001_1966 dtop001_1978 dtop001_1989 dtop001_2001
        dtop0001_1966 dtop0001_1978 dtop0001_1989 dtop0001_2001
        if spec == "psz", appendto("`tabname'.tex") rstyle(tabular);
    listtex vsp if _n == 1, appendto("`tabname'.tex") rstyle(tabular);

    listtex sec2 if _n == 1, appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab  
        dtop01_1966 dtop01_1978 dtop01_1989 dtop01_2001
        dtop001_1966 dtop001_1978 dtop001_1989 dtop001_2001
        dtop0001_1966 dtop0001_1978 dtop0001_1989 dtop0001_2001
        if spec == "alt3021", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab  
        dtop01_1966 dtop01_1978 dtop01_1989 dtop01_2001
        dtop001_1966 dtop001_1978 dtop001_1989 dtop001_2001
        dtop0001_1966 dtop0001_1978 dtop0001_1989 dtop0001_2001
        if spec == "alt3022", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab  
        dtop01_1966 dtop01_1978 dtop01_1989 dtop01_2001
        dtop001_1966 dtop001_1978 dtop001_1989 dtop001_2001
        dtop0001_1966 dtop0001_1978 dtop0001_1989 dtop0001_2001
        if spec == "alt3023", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab  
        dtop01_1966 dtop01_1978 dtop01_1989 dtop01_2001
        dtop001_1966 dtop001_1978 dtop001_1989 dtop001_2001
        dtop0001_1966 dtop0001_1978 dtop0001_1989 dtop0001_2001
        if spec == "alt3024", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab  
        dtop01_1966 dtop01_1978 dtop01_1989 dtop01_2001
        dtop001_1966 dtop001_1978 dtop001_1989 dtop001_2001
        dtop0001_1966 dtop0001_1978 dtop0001_1989 dtop0001_2001
        if spec == "alt3025", appendto("`tabname'.tex") rstyle(tabular);
    listtex vsp if _n == 1, appendto("`tabname'.tex") rstyle(tabular);

    listtex title3alt if _n == 1, appendto("`tabname'.tex") rstyle(tabular);
    listtex midrule if _n==1 , appendto("`tabname'.tex") rstyle(none);
    listtex sec3 if _n == 1, appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab 
        dtop01_1966 dtop01_1982 dtop01_1989 dtop01_2001
        dtop001_1966 dtop001_1982 dtop001_1989 dtop001_2001
        dtop0001_1966 dtop0001_1982 dtop0001_1989 dtop0001_2001
        if spec == "forbes3016", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab 
        dtop01_1966 dtop01_1982 dtop01_1989 dtop01_2001
        dtop001_1966 dtop001_1982 dtop001_1989 dtop001_2001
        dtop0001_1966 dtop0001_1982 dtop0001_1989 dtop0001_2001
        if spec == "forbes3016replace", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab 
        dtop01_1966 dtop01_1982 dtop01_1989 dtop01_2001
        dtop001_1966 dtop001_1982 dtop001_1989 dtop001_2001
        dtop0001_1966 dtop0001_1982 dtop0001_1989 dtop0001_2001
        if spec == "forbes3022", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab 
        dtop01_1966 dtop01_1982 dtop01_1989 dtop01_2001
        dtop001_1966 dtop001_1982 dtop001_1989 dtop001_2001
        dtop0001_1966 dtop0001_1982 dtop0001_1989 dtop0001_2001
        if spec == "forbes3024", appendto("`tabname'.tex") rstyle(tabular);
    listtex vsp if _n == 1, appendto("`tabname'.tex") rstyle(tabular);   

    listtex title3 if _n == 1, appendto("`tabname'.tex") rstyle(tabular);
    listtex midrule if _n==1 , appendto("`tabname'.tex") rstyle(none);
    listtex sec4 if _n == 1, appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab  
        dtop01_1966 dtop01_1978 dtop01_1989 dtop01_2001
        dtop001_1966 dtop001_1978 dtop001_1989 dtop001_2001
        dtop0001_1966 dtop0001_1978 dtop0001_1989 dtop0001_2001
        if spec == "alt3018", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab  
        dtop01_1966 dtop01_1978 dtop01_1989 dtop01_2001
        dtop001_1966 dtop001_1978 dtop001_1989 dtop001_2001
        dtop0001_1966 dtop0001_1978 dtop0001_1989 dtop0001_2001
        if spec == "alt3026", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab  
        dtop01_1966 dtop01_1978 dtop01_1989 dtop01_2001
        dtop001_1966 dtop001_1978 dtop001_1989 dtop001_2001
        dtop0001_1966 dtop0001_1978 dtop0001_1989 dtop0001_2001
        if spec == "alt3027", appendto("`tabname'.tex") rstyle(tabular);
    listtex vsp if _n == 1, appendto("`tabname'.tex") rstyle(tabular);

    listtex title3alt if _n == 1, appendto("`tabname'.tex") rstyle(tabular);
    listtex midrule if _n==1 , appendto("`tabname'.tex") rstyle(none);
    listtex speclab  
        dtop01_1966 dtop01_1982 dtop01_1989 dtop01_2001
        dtop001_1966 dtop001_1982 dtop001_1989 dtop001_2001
        dtop0001_1966 dtop0001_1982 dtop0001_1989 dtop0001_2001
        if spec == "forbes3027", appendto("`tabname'.tex") rstyle(tabular);
    #delimit cr

    qui listtex botrule if _n == 1, appendto("`tabname'.tex") rstyle(none)  
    qui listtex end if _n == 1, appendto("`tabname'.tex") rstyle(none)/*%>*/
    ****************************************************************************

    ****************************************************************************
    * Make Table Fragments (Panel D)%<
    ****************************************************************************
    clear
    set obs 1

    gen tab = "\begin{tabular}{lccccccccc}"
    gen toprule = "\toprule"
    gen midrule = "\midrule"
    gen botrule = "\bottomrule"
    gen cmidrule = "\cmidrule(lr){2-4} \cmidrule(lr){5-7} \cmidrule(lr){8-10}"
    gen vsp = "\vspace{.03in}" in 1
    gen end = "\end{tabular}" in 1

    gen c = " "

    gen title1 = "\multicolumn{10}{l}{Panel D. Growth in Top Portfolio Shares 2001--2016 (p.p.)}"
    gen title2 = "& \multicolumn{3}{c}{Top 0.1\%} &   \multicolumn{3}{c}{Top 0.01\%} &   \multicolumn{3}{c}{Top 0.001\%}"
    gen title3 = "& Fix & Pthru & C-corp & Fix & Pthru & C-corp & Fix & Pthru & C-corp"

    gen sec1 = "\textbf{Benchmark Series} & & & & & & & & &"
    gen sec2 = "\textbf{Alternative Aggregates} & & & & & & & & &"
    gen sec3 = "\textbf{Forbes Augmentation} & & & & & & & & &"
    gen sec4 = "\textbf{Robustness Series} & & & & & & & & &"

    append using `suppD'

    local tabname = "supp_series_panelD"
    
    * Top matter
    listtex tab if _n == 1 using "`tabname'.tex", replace rstyle(none)
    listtex title1 if _n == 1, appendto("`tabname'.tex") rstyle(tabular)
    listtex toprule if _n == 1, appendto("`tabname'.tex") rstyle(none)
    listtex title2 if _n == 1, appendto("`tabname'.tex") rstyle(tabular)
    listtex cmidrule if _n == 1, appendto("`tabname'.tex") rstyle(none)
    listtex title3 if _n == 1, appendto("`tabname'.tex") rstyle(tabular)
    listtex midrule if _n==1 , appendto("`tabname'.tex") rstyle(none)    

    * Formatting
    foreach v of varlist d* {
        replace `v' = round(`v', 0.1)
        format `v' %3.1f
    }

    * Sections
    #delimit ;
    listtex sec1 if _n == 1, appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab 
        dtop01_sh_hwfix dtop01_sh_pthru dtop01_sh_ccorw 
        dtop001_sh_hwfix dtop001_sh_pthru dtop001_sh_ccorw
        dtop0001_sh_hwfix dtop0001_sh_pthru dtop0001_sh_ccorw
        if spec == "baseline", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab 
        dtop01_sh_hwfix dtop01_sh_pthru dtop01_sh_ccorw 
        dtop001_sh_hwfix dtop001_sh_pthru dtop001_sh_ccorw
        dtop0001_sh_hwfix dtop0001_sh_pthru dtop0001_sh_ccorw 
        if spec == "equalrets", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab 
        dtop01_sh_hwfix dtop01_sh_pthru dtop01_sh_ccorw 
        dtop001_sh_hwfix dtop001_sh_pthru dtop001_sh_ccorw
        dtop0001_sh_hwfix dtop0001_sh_pthru dtop0001_sh_ccorw 
        if spec == "revisedsz", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab 
        dtop01_sh_hwfix dtop01_sh_pthru dtop01_sh_ccorw 
        dtop001_sh_hwfix dtop001_sh_pthru dtop001_sh_ccorw
        dtop0001_sh_hwfix dtop0001_sh_pthru dtop0001_sh_ccorw 
        if spec == "psz", appendto("`tabname'.tex") rstyle(tabular);
    listtex vsp if _n == 1, appendto("`tabname'.tex") rstyle(tabular);

    listtex sec2 if _n == 1, appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab 
        dtop01_sh_hwfix dtop01_sh_pthru dtop01_sh_ccorw 
        dtop001_sh_hwfix dtop001_sh_pthru dtop001_sh_ccorw
        dtop0001_sh_hwfix dtop0001_sh_pthru dtop0001_sh_ccorw 
        if spec == "alt3021", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab 
        dtop01_sh_hwfix dtop01_sh_pthru dtop01_sh_ccorw 
        dtop001_sh_hwfix dtop001_sh_pthru dtop001_sh_ccorw
        dtop0001_sh_hwfix dtop0001_sh_pthru dtop0001_sh_ccorw 
        if spec == "alt3022", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab 
        dtop01_sh_hwfix dtop01_sh_pthru dtop01_sh_ccorw 
        dtop001_sh_hwfix dtop001_sh_pthru dtop001_sh_ccorw
        dtop0001_sh_hwfix dtop0001_sh_pthru dtop0001_sh_ccorw 
        if spec == "alt3023", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab 
        dtop01_sh_hwfix dtop01_sh_pthru dtop01_sh_ccorw 
        dtop001_sh_hwfix dtop001_sh_pthru dtop001_sh_ccorw
        dtop0001_sh_hwfix dtop0001_sh_pthru dtop0001_sh_ccorw 
        if spec == "alt3024", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab 
        dtop01_sh_hwfix dtop01_sh_pthru dtop01_sh_ccorw 
        dtop001_sh_hwfix dtop001_sh_pthru dtop001_sh_ccorw
        dtop0001_sh_hwfix dtop0001_sh_pthru dtop0001_sh_ccorw 
        if spec == "alt3025", appendto("`tabname'.tex") rstyle(tabular);
    listtex vsp if _n == 1, appendto("`tabname'.tex") rstyle(tabular);   

    listtex sec3 if _n == 1, appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab 
        dtop01_sh_hwfix dtop01_sh_pthru dtop01_sh_ccorw 
        dtop001_sh_hwfix dtop001_sh_pthru dtop001_sh_ccorw
        dtop0001_sh_hwfix dtop0001_sh_pthru dtop0001_sh_ccorw 
        if spec == "forbes3016", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab 
        dtop01_sh_hwfix dtop01_sh_pthru dtop01_sh_ccorw 
        dtop001_sh_hwfix dtop001_sh_pthru dtop001_sh_ccorw
        dtop0001_sh_hwfix dtop0001_sh_pthru dtop0001_sh_ccorw 
        if spec == "forbes3016replace", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab 
        dtop01_sh_hwfix dtop01_sh_pthru dtop01_sh_ccorw 
        dtop001_sh_hwfix dtop001_sh_pthru dtop001_sh_ccorw
        dtop0001_sh_hwfix dtop0001_sh_pthru dtop0001_sh_ccorw 
        if spec == "forbes3022", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab 
        dtop01_sh_hwfix dtop01_sh_pthru dtop01_sh_ccorw 
        dtop001_sh_hwfix dtop001_sh_pthru dtop001_sh_ccorw
        dtop0001_sh_hwfix dtop0001_sh_pthru dtop0001_sh_ccorw 
        if spec == "forbes3024", appendto("`tabname'.tex") rstyle(tabular);
    listtex vsp if _n == 1, appendto("`tabname'.tex") rstyle(tabular);   

    listtex sec4 if _n == 1, appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab 
        dtop01_sh_hwfix dtop01_sh_pthru dtop01_sh_ccorw 
        dtop001_sh_hwfix dtop001_sh_pthru dtop001_sh_ccorw
        dtop0001_sh_hwfix dtop0001_sh_pthru dtop0001_sh_ccorw 
        if spec == "alt3018", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab 
        dtop01_sh_hwfix dtop01_sh_pthru dtop01_sh_ccorw 
        dtop001_sh_hwfix dtop001_sh_pthru dtop001_sh_ccorw
        dtop0001_sh_hwfix dtop0001_sh_pthru dtop0001_sh_ccorw 
        if spec == "alt3026", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab 
        dtop01_sh_hwfix dtop01_sh_pthru dtop01_sh_ccorw 
        dtop001_sh_hwfix dtop001_sh_pthru dtop001_sh_ccorw
        dtop0001_sh_hwfix dtop0001_sh_pthru dtop0001_sh_ccorw 
        if spec == "alt3027", appendto("`tabname'.tex") rstyle(tabular);
    listtex speclab 
        dtop01_sh_hwfix dtop01_sh_pthru dtop01_sh_ccorw 
        dtop001_sh_hwfix dtop001_sh_pthru dtop001_sh_ccorw
        dtop0001_sh_hwfix dtop0001_sh_pthru dtop0001_sh_ccorw 
        if spec == "forbes3027", appendto("`tabname'.tex") rstyle(tabular);
    #delimit cr

    qui listtex botrule if _n == 1, appendto("`tabname'.tex") rstyle(none)  
    qui listtex end if _n == 1, appendto("`tabname'.tex") rstyle(none)/*%>*/
    ****************************************************************************

end

/*%>*/
*******************************************************************************

capture program drop numbers_draft/*%<*/
program define numbers_draft

    * For facts that aren't readily available in graphs

    * Private business facts
    load_analysis_data scfrevision
    keep if year == 2016
    tabstat hwbus privccorw networth_pref [fw=wgt1B], s(sum)

    es_rank_scf, rankvar(networth_pref) othersplitvars("ccorw")
    gen top001 = es_rank >= .9999 
    tabstat ccorw [fw=wgt1B], s(sum) by(top001)

    * Social security facts
    import excel using $inputs/ssw_szz.xls, ///
        sheet("nwdbsort") cellrange(A1:C55) firstrow clear

    gen consolidate_group = cond(wealth_grpy < 6, 1, 2)
    collapse (sum) ssw_shv = exp_ssw, by(year consolidate_group)

    * Auto dealers
    load_analysis_data hetero_vals_v3
    list count profits_S_soi sale_S_soi assets_S_soi ebitd ebitd_syzz ///
    ebitd_syzz_hyb eqmultiple_sale eqmultiple_assets eqmultiple_ebitd if year == 2016 & naics_4d == 4411

end/*%>*/

*******************************************************************************
* Main
*******************************************************************************
capture program drop main
program define main

    *************************************************************
    * Set up
    *************************************************************
    !mkdir "${outputdir}/draft_final"
    log_start, f("draft-final")
    local outdir $outputdir/draft_final
    cd "`outdir'"
    *************************************************************
    * Sections
    *************************************************************
    set trace on
    set tracedepth 2

    * Creates some environment variables for preferred specs
    set_env_ds

    *********************************************************
    * Introduction
    *********************************************************
    graph_comparison_overtime

    *********************************************************
    * Hetero returns
    *********************************************************
    * Fixed
	graph_scf_portfolio_fixed
    graph_interest_participation
    graph_interest_composition
    graph_bankshares_overtime 
    graph_asset_returns
    graph_hetero_fixed
    graph_rates_overtime_fixed
    graph_return_ratio_fixed
    graph_cmd
    graph_factors_overtime_fixed
    graph_comparison_trends
    graph_mse

    * Pass-through business
    graph_aggregate_passthru
    graph_indy_passthru
    graph_passthru_hockey

    *********************************************************
    * Results and comparisons
    *********************************************************
    table_topwealth 2016
    table_portfolioshares

    graph_p9099_wealth
    graph_topwealth_cross

    graph_components 
    graph_components_overtime
    graph_uncertainty
    graph_imputation

    table_supplemental

    *********************************************************
    * Numbers
    *********************************************************
    numbers_draft

    cd $localcodedir
    log close

end

